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Section 1 

INTRODUCTION AND SUMMARY 


1.1 OBJECTIVE 

The objective of the work reported here is to develop a computer-implemented numerical 
method for predicting the flow characteristics and performance of three-dimensional jet 
exhaust nozzles. The immediate objective is to develop a method for an isolated nozzle. 
Both the internal exhaust gas flow and the external air flow with which it interacts are 
addressed. High Reynolds number turbulent flows are of primary interest for any 
combination of subsonic, transonic, and supersonic flow conditions inside or outside 
the nozzle. The method is designed to handle a variety of geometrically complex 
nozzle configurations, including nozzles that possess one plane of symmetry, two 
mutually perpendicular planes of symmetry, and nozzles with internal wedge plugs 
and external side plates. Of specific interest is the nonaxisymmetric nozzle configura- 
tion, illustrated in Fig. 1-1. The internal cross-sectional shape is circular at the 
entrance, is super-elliptical through a transition section, and becomes rectangular at 
the exit. 

The author is indebted to Dr. Dieter K. Lezius, who provided the material on turbulence 
models that appears in Section 6, and to Mr. Jacques F. Middlecoff and Mrs. Karen 
L. Neier, who performed the bulk of the effort in programming the algorithm for com- 
puter solution. 


1.2 APPROACH 

The present approach has three primary features: 

(1) A unified solution to the time-dependent Navier-Stokes equations in all 
regions of the flowfield and at all Mach numbers 

(2) A general boundary-conforming curvilinear coordinate system and 
computational grid 

(3) A fully implicit numerical method that is unconditionally stable with respect to 
the time step 
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Fig. 1-1 Sketch of Nozzle Geometry and Computational Boundaries 
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The nozzle flows in question are steady. The approach adopted here yields the desired 
steady flow solution as the time -asymptotic limit of an unsteady flow. The implicit 
numerical method is applicable to the complete three-dimensional Navier-Stokes equa- 
tions. However, in the interest of computational efficiency, we have elected to employ 
the "spatially parabolic" approximation in which the viscous transport processes 
in the general streamwise direction are neglected. The resulting equations retain all 
terms that are significant in wall boundary layers, whether separated or attached, 
and in free shear layers. They also retain all inviscid terms, so that predominantly 
inviscid regions of the flow are represented accurately. 

The ability of the approximate equations to represent wall boundary layers and free 
shear layers accurately is of great importance for predicting the nozzle performance. 
Both the nozzle discharge coefficient and the thrust efficiency are strongly affected by 
the boundary layers; the former through blockage of the mean flow as a result of 
boundary layer displacement effects (viscous-inviscid interaction), and the latter by 
skin friction drag on all internal and external walls. For subsonic nozzle flow, the 
shear layer between the exhaust jet and the external stream may also affect the nozzle 
performance. 

With this basic approach of solving the Navier-Stokes equations throughout the geo- 
metrically complex nozzle flow field, the use of some kind of implicit numerical method 
is essential for rational computation times with reasonable boundary -layer resolution. 
Virtually all nozzle internal and external flow conditions of interest are at high Reynolds 
numbers where the flow is turbulent and the boundary and shear layers are thin compared 
to a characteristic nozzle dimension. If one employs a finite-difference grid fine enough 
to accurately resolve the thin viscous layers, high computational efficiency can be 
achieved only with an implicit, unconditionally stable, numerical method. Explicit 
methods are inadequate because their conditional stability imposes severe limitations 
on the allowable time stepsize. 


1-3 



1.3 OVERVIEW OF THE REPORT 


This report presents the theoretical foundation and formulation of the numerical method. 
A manual describing the computer program, its structure, and its application comprises 
a supplement to this report. 

The basic equations and boundary conditions are presented in Section 2. The strong 
conservation-law form of the complete three-dimensional Navier-Stokes equations 
initially is formulated in Cartesian coordinates. However, a Cartesian coordinate 
system is entirely inappropriate for the geometrically complex nozzle configurations 
of interest. The equations are transformed to the equivalent strong conservation law 
form in a general system of curvilinear coordinates. One of the coordinates is oriented 
in the general streamwise direction, and the viscous transport processes in that direc- 
tion are neglected (parabolic approximation). 

The boundaries of the flow region in which the equations are solved are indicated 
schematically in Fig. 1-1. These include the nozzle walls, the symmetry planes, an 
upstream inflow boundary, a downstream outflow boundary, and one or more lateral 
outer boundaries. It is well known that the boundary conditions exert a dominant 
influence on the solution to any system of partial differential equations. The boundary 
conditions that are applied at the various boundary surfaces also are discussed in 
detail in Section 2, 

In Section 3, the physical flow region depicted in Fig, 1-1 is mapped onto three- 
dimensional rectangular computational domain by a boundary-conforming curvilinear 
coordinate transformation. The transformation, i.e., the curvilinear coordinate 
system and computational grid, is generated numerically as the solution to an elliptic 
system of partial differential equations (Refs. 4, 10, 11, and 12). The boundary 
conditions for the elliptic system are the desired locations of the grid points that lie 
on the boundaries of the physical flow region. However, the original method has the 
shortcoming that, in practice, it is difficult to control the spacing between grid points 
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in the interior of the region. It is especially important to have good control over the 
interior grid spacing in the high Reynolds number flows of interest here, where 
adequate resolution of the steep flow variable gradients within wall boundary layers 
and free shear layers demands a locally refined grid. 

The described shortcoming of the original method has been overcome. The modified 
method presented in Section 3 yields an elliptic system whose solution has the 
remarkable property that the grid point distribution in the interior of the flow region 
is controlled entirely by the selection of the grid point distribution along the boundaries. 

The numerical method employed to solve the flow equations is discussed in detail in 
Section 4. The basic method employed at interior points of the grid is an implicit, non- 
iterative alternating-direction technique (Refs. 13 and 14). The accuracy of the numerical 
solution in the physical domain depends to a great extent on the accuracy which which 
the boundary conditions are represented numerically at the boundaries of the computational 
domain. In Section 4, we formulate implicit difference equations for the boundary grid 
points that are similar to those for interior points, but that embody the boundary condi- 
tions given in Section 2. The boundary point equations have the same time accuracy as 
the interior point equations, and are compatible with the latter in their spatial order 
of accuracy. 

Section 5 presents the results of numerical experiments that have been conducted to 
test various aspects of the general method. A major portion of the experiments have 
been performed to develop and verify methods of computing the boundary conditions 
that are valid throughout the subsonic, transonic, and supersonic flow regimes. These 
experiments played a major role in the formulation of several of the boundary conditions 
given in Section 2, and of the boundary point difference equations given in Section 4. 

Turbulence models that will be used to compute the effective viscosity coefficent p 
and Prandtl number Pr in various subregions of the nozzle internal and external 
flowfield are discussed in Section 6. 
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Section 2 

EQUATIONS AND BOUNDARY CONDITIONS 
2. 1 EQUATIONS IN CARTESIAN COORDINATES 

The strong conservation-law form of the Navier-Stokes equations in Cartesian coordinates 
can be written in the compact vector notation 


K + 

q* = 

r = 

if = 

r = 

E = 


f + g + h — cr + 0 + d 

X °y z x y z 


T 

(p, pu, pv, pw, pE) 


2 T 

(pu , p + pu , puv , puw , puH) 


2 T 

(pv , pvu , p + pv , pv w , pvH) 


2 T 

(pw , pwu , pwv , p + pw , pwH) 


e + ( u 2 + v 2 + w 2 )/2 


H = E + p/p (2.1) 

where u , v , w are the velocity components in the coordinate directions x , y , z ; 
p is the density p the pressure, e the specific internal energy; and <T, if, oT 
represent the viscous stress and work terms for each coordinate direction. Presently 
available algebraic turbulence models generally assume that the first and second 
viscosity coefficients are related according to the Stokes hypothesis X = -2 p/3 . The 
viscous stress, work, and heat conduction terms then can be written in terms of the 
turbulent viscosity p and thermal conductivity k as follows. 
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cr = fJ. 


(2/3) (2u - V - w ) 

x y z 

v + u 
x y 

w + u 
X z 

I (2/3) u(2u - v - w)+v(v + u ) + w (w +u) + (k/p) T J 
*■ x y z x y ' x z' ' /rv x 


(2.2a) 


0 = p 


u + v 

y x 

(2/3) (2v - u - w ) 

y x z' 

W + V 

y z 

u (u + v ) + (2/3) v (2v 

■ y x ; y 


u - w ) + w ( w + v) + (k/M) T J 
x z y z y 


(2.2b) 


OJ = /u 


u + w 

Z X 

V + w 

z y 


(2/3) (2w z - u x - v ) 

[u (u + w ) + v (v + w ) + (2/3) w (2w - u - v ) + (k/M) T J 

^ A Ay A X y Z 

(2.2c) 

Equations (2. 1) are supplemented by the equation of state which, for a thermally and 
calorically perfect gas, may be written in the form 


p/p 


= RT = (y - 1) e 


2 / 

= c /y 


(2.3) 


where R is the gas constant, y is the specific heat ratio, and c is the speed of 
sound. 
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2. 1„ 1 Non-Dimensionalization 


The foregoing equations can be cast in dimensionless form by introducing the normalized 
flow variables 

p = p/Poo » p = p/Poo > T = T/T to , p = p/p " 

(u, v, w) = (u, v, w)/c B , e = pE/p^c^ 

(x, y, z) = (x, y, z)/L , t = c^t/L (2.4) 

where L is a reference length scale and the subscript 00 denotes some convenient 
reference state. We define the Reynolds and Mach numbers at the reference state as 

Re = p V L/p M = V /c (2.5) 

f~ 2 2 2 

where = V u » + v oo + w ® . Equations (2. 1) to (2. 3) then assume the dimensionless 
form given below, where we have suppressed the overbars for simplicity. Eqs. (2. 6) 
and all equations throughout the remainder of the report are written in dimensionless 
variables. 


% 


+ f + g + h 
x y z 


Re 1 (<T + T + oT ) 
v x y z' 


(2.6a) 


2 

c 


= p/p 


= T 


y (y - i) (e/p - v 2 /2) 


(2.6b) 
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where 


Re 

Pr 

q 

r 

r 

r 


Re /M 

OO oo 

p c A 

p 

T 

(p, pu, pv, pw, e) 

2 T 

[pu, (p/y + pu ) , puv, puw , u (p/y + e)] 

2 T 

[pv, pvu, (p/y + pv ) , pvw , V (p/y + e)] 

2 T 

[pw, pwu, pwv, (p/y + pw ), w (p/y + e)] 


0 


a 


(2/3) (2 u - v - w ) 
x y z 


P 


V + u 

X y 

w + u 
X z 

_ (2/3) u (2 u 

X 


v - w ) + v (v +u)+w(w +u)+[(y- 

y z x y x z 


ro 


0 = 




u + v 


y x 

(2/3) (2 v y 


u 

X 


W Z } 


I W + V 

y z 

Lu (u + v ) + (2/3) v (2 v 
y x' y 


u - w ) + w (w + v ) + [(y 
x z y z 


(2.7) 


1) Pr]" 1 T J 

X 


- l)Pr ] _1 T-l 

y 
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0 


u + w 

Z X 


OJ 




V + w 

z y 


(2/3) (2 w z - u x - v ) 

.u (u + w ) + v (v + w ) + (2/3) w (2 w - u - v ) + [(y - l)Pr] _1 T - 

i A 4i y ZiAV Z 


2.2 TRANSFORMATION TO GENERAL CURVILINEAR COORDINATES 


The task of solving Eqs. (2.6) numerically within a given flow region R bounded by 
a closed surface S is simplified greatly by transforming the equations to a new 
boundary-conforming curvilinear coordinate system (Refs. 1 to 5). Here, a boundary- 
conforming coordinate system is one in which the boundary surface S is a coordinate 
surface. 

We introduce a time -dependent invertible mapping transformation of R onto boundary 
conforming curvilinear coordinates (Refs. 3 to 6 ). Here the precise meaning of the 


x, y, z 


T 



£ > f] * £ » T 


£ = 4 (x, y, z, t) 

rj = tj (x, y, z, t) 

£ = £ (x, y, z, t) 

T = t 


( 2 . 8 ) 


term "boundary-conforming" is that, in these coordinates, the boundary S is composed 
only of segments of coordinates surfaces | = const. , i? = const. , £ = const. That is, 
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the surface S has the time-invariant functional representation f (£, 77 , £) = 0 
independently of time, r. The mapping otherwise may be arbitrarily defined if it is 
invertible and satisfies this constraint. 

In applying this transformation to Eq. (2. 6 a) one usually retains the dependent variables 
associated with the original Cartesian coordinate system. 

Upon transforming to £, 17 , t coordinates with the aid of the chain rule for partial 
derivatives, Eq. (2. 6 a) becomes (Refs. 3 to 7)) 


«r + r ( + g n + 


A — 1 /a A a \ 

h S ’ Re h + % + V 


q 

A 

f 


= J 


q + £ f + £ g + £ h 
t x y z 


g = v t q + V x f + V y g + V z h 

A A A ► A ► A ► 

h = L q + i f + f g + £ h 

t M x y b z 


cr = 


£ (7 + I 0 + £ cu 

s x ^y s z 


1 ? X or + r) 0 + T, z o> 


w = £ a + £ 0 + £ w 
x y z 


(2.9) 


where J is the Jacobian of the inverse transformation T 


-1 


j = d ( x , y, z ) 
a (4, V, £) 


( 2 . 10 ) 
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and where 


h - J «, • - J «x - 

A* 

£ = J £ , etc. 

y y 

(2.11) 

We observe that Eq. (2. 6a) transforms into Eq. (2.9) only if the mapping functions 
(2. 8) satisfy identically the following relations (Ref. 6) 

A A A 

£ „ + n + t 

5 x| 'xrj fe x? 

= 0 

(2.12a) 

+ V? + 

= 0 

(2. 12b) 

A A A 

+ ^ZT] 

= 0 

(2.12c) 

A A 

J T + «ti + % + 

O 

II 

< Kt* 

(2.13) 


One can verify by direct calculation that these relations hold analytically as long as 
the mapping functions possess continuous second partial derivatives. 

The mapping functions (2. 8) need not be known analytically in order to obtain numeri- 
cal solutions to the transformed Navier-Stokes equations (2.9). Once the £, i), f 
computational domain is covered by a finite -difference grid, the inverse transforma- 
tion T 1 can be defined numerically in terms of the Cartesian coordinates x(r), y(T), 
z(t) of the grid points, from which the partial derivatives x T , y T , x^, y^, etc., are 
computed by using difference formulas. The coefficients that appear in Eq. (2.9) 
then can be evaluated from the following identities (Ref. 6): 
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A 


z £ ' 

- y £ 

z 

V 

A 

= ^ 







A 

5 y 

= z 

V 

V 

- z £ 

X 

V 

A 



x £ 





A 

^Z 

= X 

V 

y r 

' X £ 

y v 

A 

^Z 

= 


y f 





A 

f x 


z - 
V 

- y T) 

z l 

A 

*t 


X + 
T 

? y y r 

+ 

A 

^Z 

z ) 

T 7 


A 

f y 

= Z ^ 

x - 
V 

- Z 

V 


*»t 

= “ ( 7) 

V l X 

X + 
T 

n y 

ly J T 

+ 

A 

77 

V 


A 

f z 

= X 4 

y„ * 

V 

- X 

V 

y ! 

A 

‘t 

A 

= - (f 

X 

X + 
T 

A 

f y y T 

+ 

A 

t 

Z 

Z T> 

(2.14) 


Care must be exercised, however, to ensure that the finite -difference representation 

A A 

of the spatial derivatives x^, y^, x^, f^, g^, etc. is such that the identities (2.12) 
hold true numerically (Ref. 7 ). Similarly, for time -dependent transformations 
associated with moving boundaries, special measures must be taken to ensure that 
Eq. (2. 13), the "Geometric Conservation Law," also is satisfied numerically (Ref. 6 ). 
The reason lies in the fact that violation of any of Eqs. (2. 12) or (2. 13) induces first- 
order numerical errors in conservation of mass, momentum and energy. 

2.3 PARABOLIC APPROXIMATION 

For the nozzle flow problems of present interest, we shall assume that both the 
Cartesian and the curvilinear coordinate systems are oriented with the x and £ 
axes running in the general streamwise flow direction. Because the dominant viscous 
effects are associated with the cross-stream directions rj, J, we may employ the 
so-called "parabolic approximation. " In this approximation, we simply neglect all 
viscous terms in the R. H. S. of Eq. (2. 9) that involve derivatives with respect to £ . 
Before doing so, however, we shall examine the structure of the viscous terms. 

A 

Consider, for example, the quantity 9 defined in (2. 9). Upon inserting a, 9, co 
as defined below Eq. (2. 7) and transforming to |, ??, £ coordinates with the chain 
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rule, one finds that the resulting equation can be split naturally into three terms, 
each of which contains derivatives in a single coordinate direction 


e = + qW + q(£) 


(2. 15) 


0 

(Vtj • V|)u + (| Vtj - 2/3 tj V|) . V, 
s x x ? 

IP \ 

6 (Vtj • V|) v^ + (| Vtj - 2/3 tj V| ) • V 

(Vtj . V|) w + (^Vrj - 2/3 ^ V|) • V 

.(Vtj • V|) [T/(y - l)Pr + V 2 /^] + [(V • V|)Vtj - 2/3 (V * Vtj)V|] • - 

(2. 16a) 


0 

(Vtj • Vtj) u + 1/3 tj (Vtj • V ) 

'/ X 'I 

9 (V) =piJ (Vtj • Vtj) v^ + 1/3 Tj y (Vtj . V ) 

(Vtj • Vtj) + 1/3 t? z (Vtj • V^) 

. (Vtj • Vtj) [T/(y- l)Pr + V 2 /2] + 1/3 (V • Vtj) (Vtj • V ). 

(2. 16b) 

0 

(Vtj • V£) u + (£ Vtj - 2/3 tj V£) • V. 

S X X s 

= (XJ (Vtj • VS) v + (£ Vtj - 2/3 tj VS) • V ? 

^ y y ^ 

(Vtj • VS) w. + (S Vtj - 2/3 tj V£ ) • v\- 

b Z Z o 

. (Vtj • VC) [T/(y - l)Pr + V 2 /2] ? + [ (V • VS )Vtj - 2/3 (V • Vtj) VS ] • V* • 

(2. 16c) 
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where standard Cartesian vector notation is used: 

V = (u.v.w) , V = I VI , V = (%>%>%) 

vn = (\’ r l y ’ , n z ) » etc. 

Equations (2.15) and (2.16) have useful symmetry properties. For example, Eq. (2.16b) 
defining 0^ can be deduced directly from «<«> in (2. 16a) by making the simple substi- 
tution £ -» tj. Similarly, 0^ in Eq. (2 . 16c) can be obtained from (2. 16a) by making 
the substitution £ -» 

The transformed viscous terms a and to each can be split in analogous fashion: 


= &(£) +£(*)+£(£) 
a = a(0 + £0 ?) + £(f) 


(2.17) 

(2.18) 


The individual terms in (2.17) can be deduced from Eq. (2.16) by making the substitu- 
tions 0 — » co, 77 -* ?. Similarly, the terms in (2.18) can be deduced from (2.16) by 
substitution of 0 -*■ o , 17 -*• £ . 

With Eqs. (2.15) through (2.18) in hand, we obtain the parabolic approximation simply 
by neglecting all £ -differentiated terms 

&£ = = 0 (2.19) 

The Navier -Stokes equations (2.9) then reduce to 
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2.4 CURVILINEAR COORDINATE GEOMETRY AND BOUNDARY CONDITIONS 


2.4.1 Coordinate Geometry 

Equation (2.20) is to be solved numerically, subject to certain admissible boundary 
conditions that must be satisfied on the various computational boundaries illustrated 
in Fig. 1 - 1 . The curvilinear coordinate system (£ , i? ,? ) is defined such that the 
computational boundaries coincide with coordinate surfaces. We shall restrict atten- 
tion to nozzle configurations that have either one plane of symmetry y = 0 , or two 
mutually perpendicular symmetry planes y = 0 and z = 0 . The two types of con- 
figuration are illustrated in Fig. 2 - 1 . The curvilinear coordinate system, which may 
be nonorthogonal , is oriented so that the £ coordinate runs in the general streamwise 
direction. For convenience, the family of coordinate surfaces £ = const, generally 
will be taken as planes normal to the Cartesian X-axis. 

The upstream inflow and downstream outflow boundaries are defined to be members of 
this family £ = const. (Fig. 1-1). The boundary-conforming cross-stream coor- 
dinates tj , f are defined so that all nozzle walls either are members of one or the 
other of the two families of coordinate surfaces r> = const, or £ = const, or are 
composed of intersecting segments of both families 17 = const, and f = const. The 
lateral outer computational boundary is defined in similar fashion as a surface 
tj = const, or f = const. , or as a composite of the two. The symmetry planes are 
also members of one of the latter families. For example, the symmetry plane in 
Fig. 2-la is represented as a surface V - const. , whereas in Fig. 2-lb, one sym- 
metry plane is a surface 1 ? = const. , and the other is a surface f = const. The 
boundary conditions that hold at the various boundary surfaces are as follows. 

2.4.2 Wall Boundary Conditions 

At the nozzle walls , the velocity vector must vanish and either the wall temperature T 

w 

or heat flux may be specified. For the latter, we restrict attention to the 
adiabatic wall case Q = 0 . 
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u = v = w = 0 (2.21) 

and 

T = T w (2.22) 

or 

-k(n • VT) = 0 (2.23a) 

where n is the local vector normal to the wall. For a wall that coincides with the 
surface v = const. , the normal vector is it = Vtj, and Eq. (2.23a) becomes 

Vr? • VT = (Vr? • Vl) T^ + (Vi? • Vij) + (V7 • 7f) T ? = 0 

However, the term involving is neglected in the parabolic approximation, and the 
adiabatic wall boundary condition is then 

(Vr? • Vtj) + (Vi? • V?) T f =0 (2.23b) 

The corresponding boundary condition for an adiabatic wall that coincides with the 
surface f = const, (see Fig. 2-la) can be deduced from (2.23b) by making the sub- 
stitution i? -*■ f , £ -*■ 7 . 

2.4.3 Symmetry Boundary Conditions 

At a symmetry plane, the normal component of velocity is an odd function of position 
relative to the plane, whereas all other flow variables are even functions. The 
Cartesian coordinates also are either even or odd at a symmetry plane. For example, 
at the symmetry plane rj = 0 in Fig. 2-la, y and v are odd functions of rj ; 
whereas x, z, and all flow variables other than v are even 

y(-T?) - -y(r?) , v(-tf) - ~v(r}) (2.24a) 
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Fig. 2-la Configuration With One Symmetry Plane 

z 

* 



Fig. 2 -lb Configuration With Two Symmetry Planes 
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x(-rj) = x(tj) , z(-T?) - Z(T?) , u ( —17) - U(t?) 

w ( ~ 7 j) = w (tj) , p(-Tj) - p(v) > etc. 

As is well known, any odd function of 17 must vanish at 17 = 0. 


(2.24b) 


This information can be used in Eq. (2.20) to evaluate the rj -differentiated terms at 
the symmetry plane 17 = 0. One need only apply the well-known theorem which 
states that the derivative of an even function is odd, the derivative of an odd function 
is even, the product of an even and an odd function is odd, and the product of two even 
or of two odd functions is even. 

It then follows from Eq. (2.24) that, for each of the vectors g, and 0^, the 

third component is even, whereas the other four components are odd functions of V. 
The derivative of the third component thus is odd, and hence must vanish at rj = 0. 



(2.25) 


The derivatives of the other four components are even, and can be evaluated in the 
one-sided sense. For example, for the first component of g, we have 



lim 

Tj “► 0 + 


g 1 (h) - g x (0) 

n 


(2.26) 


The horizontal symmetry plane J = 0 in Fig. 2 -lb can be treated in similar fashion. 


2.4.4 Freestream Boundary Conditions 

If the external freestream conditions are selected as the non-dimensionalizing refer- 
ence conditions in Eqs. (2.4), then the following boundary conditions must.be satisfied 
in the asymptotic sense at infinite distance from the nozzle: 
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p 


1 


p = 1 


u 


(2.27a) 



It is assumed that the x-axis is aligned with the freestream velocity vector, so that 
V oo = W oo = 0. In practice, boundary conditions (2.27) must be applied at lateral 
outer computational boundary located a finite distance from the nozzle (Fig. 1-1), and 
it is not always permissible to impose the remaining freestream conditions v = w = 0 
at that boundary. Numerical experiments have shown difficulties in converging to a 
time-asymptotic steady state solution for subsonic and transonic flow when the latter 
conditions are imposed (see Section 5). In general, the velocity components v and 
w must be computed at the boundary as well as in the interior flow region. 


2.4.5 Upstream Inflow Boundary Conditions 


To discover what boundary conditions are permissible at the inflow boundary, we 
appeal to the quasi-one-dimens ional unsteady method of characteristics applied to 
Eq. (2.20). According to the method, the five -component vector equation (2.20) is 
regarded as a hyperbolic system of partial differential equations in r and £ . The 
tj and ? differentiated terms are treated simply as inhomogeneous "forcing function" 
terms. The analysis is cumbersome; we shall omit the details and merely summarize 
the major results for the case where £ = £(x) only, i.e. , where the inflow boundary 
is a plane normal to the x-axis. 


The hyperbolic system of five scalar equations has three distinct characteristics 
X. = d£/dr, one of which is of multiplicity three 



X 4 = (u + c) £ x 


(2.28) 


X 5 = < u - c) fx 
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where £ x > 0 for any non-singular coordinate transformation (2. 6). The charac- 
teristic form of the system is 


< P r 


V 

r 

w 

r 

+ X 3 P {> ' c2(p r 
+ X 4 P|) + pc(u T 

+ X 5 P^) - pc(u r 


+ X v = 
A 1 £ 

f l 

(2. 29a) 

+ x x w £ * 

f 2 

(2.29b) 

+ v £ > = 

f 3 

(2.29c) 

+ X 4V = 

f 4 

(2. 29d) 

+ v { > * 

f 5 

(2. 29e) 


where f^ , , ... are cumbersome functions of the flow variables and of the r\- and 

f- differentiated terms in (2.20), and where the operators 


_d_ 

dr 


+ X 


_d_ 
i d£ 


represent directional derivatives along the characteristics. 


If u is supersonic (u > c), all five characteristics (2.28) have positive slope in the 
£ — r plane; i.e. , all characteristics are directed downstream. We conclude that all 
flow variables may be specified a priori as inflow boundary conditions whenever u > c. 


For subsonic inflow u < c, all characteristics remain of positive slope except X 

o 

which has negative slope. Information propagates upstream along that characteristic 
from the interior of the flow region toward the inflow boundary, and the last of 
Eqs. (2.29) yields one relation among the five independent flow variables. Thus, four 
independent quantities may be specified as boundary conditions. Of these, one sees 
from Eqs. (2. 29a ,b) that two must specify the transverse velocity components or, 
equivalently, the velocity vector's direction cosines. It appears that some freedom 
exists in the choice of the remaining two boundary conditions. However, Eq. (2.29c) 
may be recognized as the isentropic relation between pressure and density perturbations 
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along a "streamline." This suggests that the entropy or, equivalently, the stagnation 
pressure, may be specified as an inflow boundary condition. On physical grounds, we 
elect to specify the stagnation temperature as the remaining inflow boundary condition. 

The boundary conditions on the stagnation pressure P 0 (£ = 0 »*?*£) and stagnation 
temperature T q (£ = 0 , rj,f) provide two independent algebraic relations among the 
local values of the flow variables at each point of the inflow boundary \ = 0. These 
relations can be represented in terms of the particle -isentropic flow equations 

T q /T = 1 + ( 7 - 1) M 2 /2 

p/p = (T/T jtAt- 1 ) 
o o 

2 2 

where M = V /T is the local Mach number written in terms of the dimensionless 
velocity and temperature [Eqs. (2.4)1. With the aid of the equation of state (2.6b), 
these relations can be written in the form 


T Q = (7 - 1) [7 e/p - (7 - l)V 2 /2] (2.30a) 

P 0 = pT o j 7 |l/(7 - 1) - e/pT o j|“ 1/(7_1) (2.30b) 


The described inflow boundary condition results for the three-dimensional "parabolized” 
Navier-Stokes equations (2.20) are essentially the same as those presented in Ref. 8 
for two-dimensional inviscid flow. 


2.4.6 Downstream Outflow Boundary 


Arguments similar to those given above for the upstream inflow boundary may be used 
at the outflow boundary. For u > c, all characteristics have positive slope and 
emanate from the interior toward the boundary. The flow variables at the boundary 
are completely determined by the interior flow, and no boundary conditions may be 
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specified. For u < c , the single characteristic emanates from the boundary, 
and one boundary condition may be specified. It is common to specify the pressure 
as equal to its freestream value (Refs. 8 and 9). 

p = 1 when u < c (2.31) 

However, we shall employ this boundary condition only when dealing with the internal 
flow in an isolated nozzle; that is, when the outflow boundary is located at the nozzle 
exit plane. For general cases involving both internal and external flow where the 
outflow boundary is located downstream of the nozzle exit, we shall specify no outflow 
boundary conditions. All flow variables at the outflow boundary will be computed from 
Eqs. (2. 20) as described in section 4. 2. 2. This approach has the following 
justification. 

The system (2.20) is not strictly hyperbolic because the viscous terms on the R.H.S. 
have parabolic qualities. Parabolic systems may be viewed as possessing degenerate 
vertical characteristics d£/dr = 0 along which information is propagated instan- 
taneously from the lateral outer boundary. Freestream pressure is imposed at the 
latter boundary [see Eq. (2.27a)l, and we rely on the parabolic behavior of the 
equations to relay the freestream pressure information to all points of the outflow 
boundary. The numerical experiments described in Section 5 attest to the validity of 
this approach. 
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Section 3 

GENERATION OF COMPUTATIONAL GRIDS 


The generation of computational grids suitable for carrying out accurate numerical 
solutions to the three-dimensional Navier-Stokes equations is currently the subject 
of intensive research. For a wide class of nozzle configurations that includes the 
configuration depicted in Fig. 1-1, a three-dimensional grid can be built up by 
constructing a sequence of two-dimensional grids in successive cross-sectional 
planes. For each cross-sectional plane between the upstream inflow boundary and 
the nozzle exit, a two-dimensional grid must be generated both for the flow region 
internal to the nozzle and for the flow region outside the nozzle. The two flow regions 
are illustrated in Fig. 3-la for a nozzle whose interior and exterior walls have a 
super-elliptical shape described by the equation 


(y/a) N + (z/b) N = 1 (3.1) 

where a,b, and N may vary with position, x, along the nozzle. For the present, 
we shall restrict attention to nozzle configurations that have two mutually perpendicular 
planes of symmetry, and shall assume that the Cartesian coordinate system is oriented 
so that the symmetry planes coincide with the coordinate planes y = 0 and z = 0 , 
respectively. 

Two substantial advantages are gained by introducing a boundary-conforming 
curvilinear coordinate transformation that maps the physical flow region in the cross- 
sectional plane (y, z) onto a rectangular computational domain (tj.J). First, a 
uniform rectangular grid can be employed in the (ri , £) domain that simplifies the 
numerical finite-difference representation of the flow equations. Second, the boundary 
conditions can be computed more easily and more accurately. Because of the boundary- 
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conforming character of the mapping transformation peripheral grid points lie directly 
on the boundaries. This avoids the tedious and inaccurate interpolation or extrapolation 
procedures that would be required to compute or to impose boundary conditions when 
the boundaries lie between grid points. 

Thompson, Thames, and Mastin (Ref. 4) have given a general method for generating 
such boundary-conforming coordinate transformations. In the method, one computes 
the (y , z) coordinates of the grid points in physical space as solutions to an elliptic 
system of partial differential equations. The equations are formulated in a rectangular 
computational domain (p , f ) , and are solved numerically over a uniform, rectangular 
grid in that domain. 

Figure 3-2 shows an example of such a mapping of the region interior to the nozzle of 
Fig. 3-1 onto a rectangular (rj , £) computational domain. Points A, B, C of Fig. 3-1 
map onto the corners A' , B', C’ of the computational domain in Fig. 3-2, where the 
point B is placed at the location where the curvature of the nozzle wall is greatest. 

This mapping produces a quasi-rectangular grid in the physical domain of Fig. 3-1 
when the Thompson, Thames, Mastin (TTM) method is applied to generate the grid 
numerically. The TTM method employs the following inhomogeneous Laplace equations 
as the generating system 


V + v * 

E yy + l zz = q< "' l> (3 ’ 2) 

These equations are transformed to f , r? coordinates by interchanging the roles of 
dependent and independent variables. This yields an elliptic system of quasilinear 
equations 

ay w - 2 %£ + yy e c ’ - j2 (Py „ + Qy i > 

az m - 2 % t + y ht ’ _j2( % +< V 
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where J denotes the Jacobian of the transformation 


(3.3) 


J 


9 (y. z ) 

9 (V, £) 


V £ 


y z 
y £ v 


(3.4) 


Equations (3.3) are solved numerically on a uniform, rectangular grid Ar] , At. 

For a simply connected physical domain such as the nozzle interior, Dirichlet boundary 
conditions may be specified over the entire closed boundary O' A' B' C' O' of the 
computational domain. The boundary values are the (y,z) coordinates of the physical 
grid points that correspond to each of the mesh points (rj^ , £^) on the boundary of the 
rectangular computational domain. 

The physical grid points may be spaced as desired along the boundaries of the flow 
region. However, the major shortcoming of the TTM method is that, in practice, it 
is difficult to control the spacing between grid points in the interior of the flow region. 
The interior grid spacing is governed primarily by the elliptic system (3. 3) itself. 

The boundary values have a strong influence on the grid spacing only in the immediate 
neighborhood of the boundaries, despite the strongly elliptic character of Eqs. (3.3). 

It is especially important to have good control over the interior grid spacing in the high 
Reynolds number flows of present interest, where adequate resolution of the steep 
velocity gradients within wall boundary layers and free shear layers demands a locally 
refined grid. 
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Various forms of the source terms P , Q in Eqs. (3. 3) have been devised that contain 
adjustable parameters and that provide some measure of control over the interior 
grid spacing. This ad hoc approach has been used with success in particular 
applications (Refs. 10, 11, and 12), but the forms of these terms and the values of the 
adjustable parameters are highly problem-dependent. No generally useful method 
has yet been devised that can be relied upon to provide an effective means of grid 
control for a wide class of boundary geometries. 

An example is given in Figs. 3-3 and 3-4 to illustrate the difficulty of controlling the 
grid in the nozzle interior. Figure 3-3 displays a quasi -elliptical grid generated by 
the TTM method by using simple Laplace equations (P = Q = 0). The parameters of 
the super-ellipse arc CD are N = 1.55, b/a = 1.4. The quasi-elliptical shape of the 
grid was achieved by specifying the boundary values so that the regular point B on 
the lower boundary AC of the physical domain maps onto the lower right-hand comer 
point of the rectangular computational domain of Fig. 3-2. The point B plays a role 
akin to that of the focus in the classical system of orthogonal elliptical coordinates. 
Equally spaced boundary values (y,z) were specified along each of the boundary 
segments AB, BC, CD, and DA. The resulting grid is smooth and regular in the 
interior, and would be suitable for inviscid flow computations. However, this grid 
would be incapable of resolving the nozzle wall boundary layer in a viscous flow. In 
an attempt to obtain a locally refined grid near the nozzle wall, an exponential dis- 
tribution of boundary values was specified along the boundary segments AB and OC. 
The resulting grid, displayed in Fig. 3-4, is wholly unsatisfactory. Evidently, the 
influence of the boundary values does not penetrate very deeply into the interior, 
in spite of the elliptic character of Eqs. (3. 3). The interior grid is affected primarily 
by the generating equations, rather than by the boundary values. 

Equally poor results are obtained whenever highly stretched boundary values are used 
with P = Q = 0. A further example for a quasi-rectangular mapping of a high-degree 
super-ellipse (N = 10) is displayed in Fig. 3-5. 


3-5 



SUPER ELLIPSE UNCLUSTERED GRID 
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Fig. 3-3 Quasi-Elliptical Interior Grid Generated by Equally Spaced 
Boundary Values and P = Q = 0 . N = 1. 55 , b/a = 1.4 


SUPER ELLIPSE 


Z 



0 1.4 


Fig. 3-4 Quasi-Elliptical Grid Generated by Exponentially Spaced Boundary Values 
on Segments AB, OC With P = Q = 0 


3-7 


SUPER ELLIPSE 



Fig. 3-5 Quasi-Rectangular Grid Generated by Exponentially Spaced Boundary 
Values and P = Q = 0. N = 10, b/a = 1.4 
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The described shortcoming of the TTM method has been overcome by devising general 
source terms P and Q that are computed from the Dirichlet boundary values, and 
that vary continuously throughout the computational domain. This results in an elliptic 
generating system whose solutions have the remarkable property that the grid point 
distribution in the interior of the physical flow region is controlled entirely by the 
user's selection of the grid point distribution along the boundaries. This is accomplished 
by choosing source terms such that Eqs. (3.3) possess exponential solutions, although 
the terms themselves are not exponentials functions. The terms contain arbitrary 
parameters that are evaluated locally at the boundaries of the computational domain 
by using a finite-difference representation of the limiting form of the system (3.3) that 
is valid at the boundaries. Interpolation of these parameters into the interior from the 
boundaries then defines the parameters at each mesh point of the computational domain. 
Numerical solution of Eqs. (3. 3) by standard successive line over-relaxation (SLOR) 
then results in a grid throughout the physical domain that is controlled entirely by the 
selection of grid points on the boundaries of that domain. 

The source terms have the form 


P = (P 07, t) (t? 2 + J? 2 ) 

y z 

Q = 07, ?) U 2 + £ 2 ) (3.5) 

where the parameters (p , ip are yet to be specified. Upon introducing these terms, 
Eq. (3. 3) assume the form 

a + %) - 2 %£ + 7 (% + (bfc) = 0 
a ( z vv + <P Z V ) - 2 PZjjg + y (Zgg + = o (3.6) 
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One can verify easily that Eqs. (3. 6) possess exponential solutions if the parameters 
(p , ip are locally constant. Given a set of boundary values (y,z) on the boundary of 
the computational domain, one can determine the parameters by requiring that the given 
boundary values satisfy appropriate limiting forms of Eqs. (3. 6) along the boundary 
of the computational domain. For example, along either of the vertical boundaries 
tj - const, in Fig. 3-2, the limiting form of Eqs. (3.6) is 


Along 17 = const. 


Along £ = const. 


y t£ + * y £ 

= 0 

(3. 7a) 

Ht + * z £ 

- 0 

(3.7b) 

boundaries 

£ = const, is 


y-nr, + 

= 0 

(3.8a) 

?777 77 



z + </?z 

7JT) 7} 

= 0 

(3.8b) 


For the quasi -rectangular mapping of the nozzle interior shown in Fig. 3-1 onto the 
computational domain shown in Fig. 3-2, a finite-difference representation of 
Eq. (3. 7b) can be used to compute the parameter ip locally at each mesh point along 
the vertical boundaries rj = const. 



ip (v, t) ■■ 

= - z c/ z e 

(3.9a) 

1! 

(Z k, f + 1 ' 

■ \l-l >/ZAe 

(3.9b) 

fgV t 

(z k, i+i 


(3.9c) 


The parameter <£(77, f) is evaluated similarly from Eq. (3. 8b). 
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For more complicated mappings one or the other of Eqs. (3. 7) may lead to an in- 
determinate equation for ip. For example, the quasi-elliptical mapping shown in 
Fig. 3-3 is such that Eq. (3.7b) is indeterminate along the right-hand vertical boundary 
of the computational domain of Fig. 3-2, and Eq. (3. 7a) must be used in place of 
(3. 7b) to evaluate ip. For an arbitrary nonsingular mapping, one can show that one 
or the other of Eqs. (3. 7) is always determinate at each point of the boundary, and 
defines a unique value of the parameter ip at that point. This follows from the fact 
that, along the boundary, y and z are related by the equation f(y,z) = 0 that 
describes the geometric shape of the boundary curve in the physical domain. Thus, 
in principle, either of Eqs. (3. 7) can be used to compute the parameter ip . In practice, 
the best results are obtained by employing Eq. (3.7a) locally at points where y^ — , 

and Eq. (3.7b) locally at points where y < z^. A similar statement holds for the 
parameter (p in Eqs. (3.8). 

Once the parameter ip is defined at each mesh point of the vertical boundaries 
77 = const, in the computational domain, its value at interior mesh points can be 
computed by linear interpolation along horizontal mesh lines f = const. Similarly, 
ip is computed by interpolation between the boundaries f - const, at which it is defined 
by Eqs. (3. 8). Equations (3. 6) then are solved by SLOR iteration in the computational 
domain to generate the grid in the physical domain. 

Figure 3-6 shows the physical grid generated by the described technique for the interior 
and exterior flow regions of the super-elliptical nozzle geometry illustrated in Fig. 3-1. 
The grid, highly refined near the nozzle walls to resolve the wall boundary layers, was 
generated simply by inputting an exponential grid point distribution along each boundary 
of the computational domain illustrated in Fig. 3-7b. Observe the smooth, regular 
character of the grid, the near-orthogonality of the grid lines to the nozzle wall, and 
the nearly uniform spacing between the wall and nearest wall-like grid line over the 
whole length of the wall. 
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SUPER ELLIPSE TEST CASE BLOWUP OF MOST CONGESTED AREA 


1.0 4 - 



Fig. 3-6b Magnified View of Grid Near Corner of Nozzle 
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Although the grid illustrated in Fig. 3-6 was generated using an exponential distribution 
of grid points along the physical boundaries of the flow region, the present technique is 
valid in general for any arbitrary boundary point distribution. In essence, the numerical 
evaluation of the parameters (p, ip at the boundaries simply is a convenient way to 
construct a curve -fit to the boundary values with local embedded exponential functions. 
The interpolation of these parameters into the interior of the computational domain 
simply extends the range of the curve-fit parameters into the interior. The elliptic 
equation system then merely provides a reliable, automatic means for translating the 
parameters into a local exponential curve-fit at each interior point that reflects the 
boundary value distribution, and that has the properties of regularity and montonicity 
required of non-singular coordinate transformations. 

For an arbitrarily selected cross-sectional plane x * const. , the described procedure 
automatically generates a boundary-conforming curvilinear coordinate transformation 
T [Eq. (2. 8)] that maps the flow region bounded between the symmetry planes and 
the lateral outer boundaries onto a rectangle 0 ^ 77 < tj ^max * 

latter then is subdivided into a rectangular grid of uniform spacing A 77 , a£ . The 
coordinates of the grid points are then 


V, = (k - 1) At? 
k 


1 < k < K 


(3.9a) 


£« = (*- 1) AS 


1 < f < L 


(3.9b) 


Repetition of the procedure for a sequence of cross-sectional planes x^ defines a 
three-dimensional transformation that maps the entire nozzle flow region of Fig. 1-1 


onto a rectangular solid 0 ^ ^ | 


0 <77 , 0 

max ’ 


£ ^ L 


max' max' ' “max 

grid spacing A| then completes the definition of the computational space 


A uniform 


= (j - 1) A£ 1 < j < J 


(3. 9c) 
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Without loss of generality, we may choose 


M = Arj = A£ = 1 (3.10) 

so that the triplet of indices (j, k, 1) describes uniquely the coordinates of each grid 

point in the computational space. For general mappings, each computational grid 

point has an image point (x. , . , y. , . , z. , ) in the physical space illustrated 

J»k,* J»k,x J>K,* 

in Fig. 1-1. However, note that the mapping used here employs cross-sectional 
planes. Consequently, the x-coordinate of a grid point depends only on the index j . 
Furthermore, we shall assume that the physical grid is stationary; i.e. , that the 
mapping functions (2. 8) are independent of time t. Time-dependent mappings are 
useful in problems that involve free boundaries such as external shock waves (see 
Ref. 6), but are unnecessary in the present application. 


The coordinates x^ define the physical location of the cross-sectional planes, and 

may be selected arbitrarily. In each such plane j , the coordinates y. ^ ^ and 

z. , of the grid points in physical space are computed numerically from the 
J » k , t 

elliptic system (3. 3). The manner in which these numerically-generated coordinate 
transformations are used in flowfield computations is explained in Section 4, which 
presents the difference equations derived from the flow equations that have been 


summarized in Section 2. 
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Section 4 

IMPLICIT DIFFERENCE EQUATIONS AND SOLUTION TECHNIQUE 


The Navier-Stokes equations (2.20) are to be solved numerically in the computational 
space defined in Section 3, subject to the boundary conditions given by Eqs. (2.21) 
through (2.31). We shall employ an implicit numerical method to avoid the severely 
limited time stepsize that would be imposed by the numerical stability bound of an 
explicit method. For the latter type of method, one can show that the maximum stable 
time step allowed by an explicit scheme is proportional to the smallest distance 
between any two grid points in physical space, whenever that grid spacing is suffi- 
ciently fine to resolve accurately the thin viscous wall boundary layers and free shear 
layers in the nozzle flowfield. The resultingly great number of time steps needed to 
attain convergence to the desired steady-state solution easily can lead to impractically 
long computer runtime. Properly formulated implicit methods are free of such severe 
numerical stability criteria. 

The numerical methods developed by Beam and Warming (Ref. 13) and by Briley and 
MacDonald (Ref. 14) basically are similar in their use of implicit time-differencing 
and of alternating direction techniques (spatial operator factorization) for the multi- 
dimensional Navier-Stokes equations. Beam and Warming first applied the method to 
the conservation-law form of the full Navier-Stokes equations for two-dimensional 
laminar flow in Cartesian coordinates. The method was combined with the general 
mapping in Eq. (2.8) by Steger (Ref. 3) for two-dimensional plane flow; by Kutler, 
et al. (Ref. 15) for axisymmetric flow; and by Pulliam and Steger (Ref. 7); and by 
Thomas and Lombard (Ref. 6) for three-dimensional flow. However, the latter four 
references employed the "thin-layer" approximation to the Navier-Stokes equations in 
which the viscous terms associated with only one of the three coordinate directions are 
retained. The thin layer approximation is obtained by dropping the first three terms 
from the R.H.S. of Eq. (2.20). We shall apply the Beam-Warming method to Eq. (2.20) 
itself. We derive first the difference equations for interior grid points. The boundary 
conditions require special treatment, and will be dealt with later in Section 4. 2. 
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4. 1 INTERIOR GRID POINTS 


In Ref. 13, Beam and Warming approximated the spatial derivative terms in the 
Navier-Stokes equations by standard second-order accurate central difference 
operators, and employed general implicit time-differencing of either first- or 
second-order accuracy, depending upon the choice of a parameter. We shall employ 
the first-order accurate Euler-implicit time-differencing, since we seek only the 
asymptotic steady-state solution. The accuracy of the steady-state numerical 
solution is governed by the spatial difference operators (Ref. 13). 

4.1.1 Time Differencing 


The first-order accurate implicit time-differenced form of Eq. (2.20) is 


Aq + /af a| ah \ n+1 

At + y$T + dr? + 3f J 


R 


67] 9f 


n+1 


# + £.£<>» 
drj 3^ 


(4.1) 


where the superscript n denotes the time level, A is the classical forward time 
difference operator 

. * a n+1 An 
Aq = q - q 


» n+1 n 

Ar = r - r 


(4.2) 


and, for later convenience, we have switched from the subscript notation to the 
standard partial differential operator notation for spatial derivatives. 

Note that the viscous cross-derivative terms — the terms in the second set of brackets 
in Eq. (4. 1) — are evaluated explicitly at time r n , whereas all other terms are evalu- 
ated implicitly at the advanced time r n+ * = r n + Ar (Ref. 13). As shown in the latter 
reference, this treatment of the cross-derivative terms does not compromise the 
stability of the numerical method. The remaining terms are locally linearized in time 

A A 

by using local Taylor series expansions. For the inviscid flux vectors f and g, the 
linearization is 
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(4.3) 


^n+1 


*n+l 

g 


? n +1 
h 


f n 


An 

g 




r.n a 

+ F Aq + • • • 

A a 

+ G Aq + • • • 
+ H Aq i- • • * 


where F,G,H, are the Jacobian matrices 

F = dt/d q G = 3g/3q H - 3h/3q 


It follows from Eqs. (2. 9) that the Jacobian F is given by 

F = £ I + £ F' + £ G' + £ H' 
t x y z 


(4.4a) 


(4.4b) 


F' = 3f/3q G' = 3g/3q H' = 3h/3q’ (4.4c) 

I = Identity matrix ( 4 . 4d) 

The corresponding expressions that define matrices G and H can be obtained from 
(4.4b) by the substitutions £->-77 and £ -*■ J, respectively. 

The Taylor series expansions given in Eq. (4.3) are valid only if the coordinate 
transformation ( 2 . 8 ) is independent of time 

K = = h = 0 < 4 - 4e > 

This implies that the grid is stationary in physical space (x,y ,z). The reader is 
referred to Ref. 6 , where the general case is treated for moving grids associated 
with time-dependent coordinate transformations. 

The viscous terms are linearized similarly as follows: 
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n 


(4. 5a) 


e 



n+l 



—**•*> *\ 

+ R Aq + 


co 


a) 


n+1 


co 


(0 


0 n a 

S Aq + 


(4.5b) 


where R,S are matrices whose elements are differential operators. The notation 
such as RAq thus is to be interpreted not as the product of a matrix R and a vector 
Aq, but rather as a matrix operator that operates on the vector Aq* The expressions 
that define the operators R , S are derived in Appendix A. 


The expansions (4.3) and (4. 5) allow Eq. (4. 1) to be written in the operational form 



where the linear operator notation is to be interpreted as follows: 



(4.7) 


That is, the spatial differential operator acts upon the product F n Aq and not upon the 
matrix F n itself [c£. Eqs. (4. 5) and (A. 6), and Ref. 13]. TheR.H.S. of Eq. (4.6) 
consists of the spatially differentiated terms in Eq. (4. 1), evaluated at time level n 


An j 9f 9 

r = ^ 


a -l*m 
g - Re 


Re 


h - Re^ (T?) 


_L gfo> + A *n> 


3rj 




(4.8) 


The various terms in Eq. (4. 8) are grouped in a fashion that is convenient for the 
spatial differencing to be performed later. 


4-4 



The final form of the time-differencing algorithm is obtained by noting that the 
spatially three-dimensional operator on the L.H.S. of (4.6) can be factored into a 
product of one-dimensional operators, with a remainder term that is of the same 
order as the temporal truncation error of (4.6) itself (Ref. 13). The factored 
operator equation is 



(4.9) 


Although the three operators are not commutative, the truncation error is of second 
order in At for any arrangement of the three factors. 


Once the spatial derivative operators are replaced by finite-difference operators, the 
solution vector q is obtained through the ADI (Alternating Direction Implicit) 
sequence 



An 

r 



where the operators in parentheses in the second and third factors of Eq. (4. 9) have been 
rearranged in a form more convenient for spatial differencing. 


4-5 



4*1.2 Spatial Differencing 


The spatial differencing of Eqs. (4.10) can be represented conveniently in terms of 

1/2 

classical difference operators. We shall employ the half-mesh shift operator E 

/ 1/2 -l/2\ 

and two central operators: the averaging operator fj. = I E + E f/2 , and the 

1/2 - 1/2 ' 

central difference operator 5 = E - E . For a mesh function f defined at 
mesh points j , the operators are defined as follows (Ref. 16): 


f j = f j±l/2 

(4.11a) 

Mf j = (Vl/2 + f j-l/2 )/ 2 

(4.11b) 

5f j = f j+l/2 " Vl/2 

(4.11c) 


However, we must deal with three-dimensional mesh functions such as 
confusion about the spatial direction in which a given operator acts, we 
subscript to indicate that direction. For example, we have 


f., , . To avoid 

Jkl 

shall append a 


f jkl 


( f J. 


k+1/2,1 + f j , k-l/2,l)/ 2 


(4,12) 


The single subscript also will serve to distinguish between the averaging operator and 
the viscosity coefficient, both of which are represented by the symbol /i . 

We consider first the spatial differencing of the term f n in the first of Eqs. (4.10) as 
defined in Eq. (4.8). The superscript n indicating the time level will be suppressed 
for brevity. The first-order spatial derivative operator 3/9£ of the first term in 
braces in Eq. (4.8) is represented as follows by a central difference formula that is 
accurate to second order in the spatial grid spacing . 


H ~ ("j w/* s ( f , 


j+1 , k , 1 f j-!» k 


,i)/ 2A f 


(4.13a) 


A£ = 1 
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Ths increment A£ will be omitted from this and from subsequent spatial difference 
formulas by virtue of Eq. (3.10). The derivative operators 3/3 tj , 3/3J of the second 
and third terms in braces in Eq. (4.8) are represented by the corresponding central 
difference operators 


JL 

3t? 


JL 

as 


g - Re 


1 $<r>' 


h - Re 


1 £«] 


“k \ [* 

"i h l 


Re 1 

„-l A 

Re co 



jkl 


jkl 


(4. 13b) 


(4.13c) 


The same difference operators are used to represent the derivatives that enter into the 
vectors 0^ and co^. For example, the first term in the second component of the 
vector 0^ is [Eq. (2.16c)]. 


M J (Vi? * V?) Uj. 


qu J ( 7 7 f +TJ f +17 f 

1 w x J x yy ‘z J z 


jjkl 


^1 6 1 


U. M 

jkl 


(4.13d) 


One sees from Eq. (2. 16b) that the first term enclosed by the last set of brackets in 
Eq. (4.8) involves repeated derivatives with respect to rj . Botn the outer derivative 
3/3 tj of the term in question and the inner derivatives that enter into 0 ^ in 
Eq. (2.16b) are represented by the central difference operator 5. 

K 


_3_ 

3tj 



(4. 14a) 


Thus, quantities such as 8/817 (MJ IVhl 2 3U/3tj) that appear when the term dd^/drj is 
written in expanded form are differenced as follows: 


i?(“ J| ™ | 2 !f)~ s k (^ivil 2 s k u) 


(4.14b) 


The last term in Eq. (4. 8) is differenced similarly in terms of the operator <5^ 


JL 

3£ 
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(4. 14c) 



Each spatially differentiated term on the left of Eqs. (4. 10) involves either an ordinary 
matrix such as F or a matrix of spatial differential operators such as R„ Each such 
matrix is the Jacobian of one of the terms that comprise the quantity r in (4. 8). For 
consistency, the spatial derivative operators that appear in each term on the left of 
(4.10) are represented b}' the same difference operators that are used for the corre- 
sponding term of r in (4. 8). For example, in the first of Eqs. (4.10), the term 

A 

involving the Jacobian of the inviscid flux vector f in Eq. (4. 8) is differenced as 
is (4. 13a). 

>1.8. F (4.15) 

Similarly, the terms (0/Dr])G and (3/3£)H are differenced as in Eqs. (4.13b) and 
(4.13c), respectively. Finally, Eqs. (4. 14a) and (4. 14c) are used for the last term in 
parentheses in each of the second and third of Eqs. (4. 10). 


When the curvilinear coordinate system and the computational grid are generated 
numerically as described in Section 3, the grid generation procedure yields only the 
x,y ,z coordinates of each grid point in physical space. The "metric coefficients" 

, etc. then must be computed from Eqs. (2.14). This can be accomplished by 

x y 

using the central difference operator fiS to represent the spatial derivatives. The 
final difference equations for interior points then take the operational form 


[i + A Vj F r A ’iki = % 


(4. 16a) 


I + At - Re 


k~k 


V)]" A 5* k l - A ^! 


(4.16b) 


[■ 


+ A r - Re S^S^ 

A ^ki 

- Ski 

(4. 16c) 


a n+1 
q jkl 

= Sjh + A V 

(4. 16d) 
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where 



~ At |mj 5^ f + M k § k [g - Re 1 


ij 5j h - Re 1 gj^ - Re 1 « 


5 , 0 

k 


0 ?) 


+ 5. 




(4.17) 


and where it is understood that the spatial derivatives involved in 0^ and co^ are 
differenced as described following Eq. (4. 13c), and those involved in 6 ^\ R , 

and S are differenced according to Eq. (4.14). 


Each step of the ADI sequence (4. 16) involves the solution of a linear system of equa- 
tions having a block-tridiagonal coefficient matrix. For example, upon applying the 
operator on the L.H.S. of Eq. (4.16a), we obtain 


4 r ( f " a ^-*i) + + % • 


2 <; j s J - 1 


(4.18) 


where the subscripts k,l have been suppressed. The block-tridiagonal structure is 
readily apparent, inasmuch as the Jacobian matrices F? + ^ (the blocks) are 5X5 
square matrices. 


However, the system of Eq. (4. 18) is incomplete; it involves only J-2 equations among 
the J unknown values of Aq** . The missing equations are those for the boundary points 
j = 1 , J . To close the system, Steger (Ref. 3) and Pulliam and Steger (Ref. 7) simply 
assume Aq|* = Aqjf* = 0 ; i. e. , that the flow conditions at the boundaries do not 
change over a time step. After the sequence (4.16) is completed for interior points, 

— V 

they then update the values of q at the boundaries by ad hoc methods that are intended 
to satisfy the boundary conditions at steady state. This time-lagging approach implies 
only a zero-order time accuracy at the boundaries. The approach does not necessarily 
degrade the accuracy of the steady-state solution, depending on the updating procedures, 
but it may retard the rate of convergence to steady state. The numerical experiments 
described in Section 5 suggest that this is the case. 
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Rather than use the time-lagging approach, we follow the approach taken in Ref. 15, 
and formulate implicit difference equations for the boundary grid points that are 
similar to those for interior points, but that embody the boundary conditions of 
Section 2.4. 

4. 2 BOUNDARY POINTS 

In this subsection, we formulate boundary point-difference equations that have the same 
time accuracy as the interior point equations, and that are compatible with the latter 
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Fig. 4-la Sketch of Computational Grid in a Plane t = Constant 


Z 




at a cell face is defined as the average of 


In the finite volume method, the value of f 
values at the two neighboring grid points 


f j + l/2,k,f (f jkj2 + f j + l,M )/2 M j f j + 1/2,M 


(4.20b) 


When the left-hand equality of (4.20b) is substituted into that of (4.20a), the resulting 
equation is identical to Eq. (4. 13a); namely, the central difference that is obtained by 
applying the operators 6^ and / j .. in sequence. Similar results hold for all the spatially 
differentiated terms of Eq. (4. 1). Furthermore, by the mean value theorem, the first 
term in the volume integrated version of (4. 1) represents the time derivative of the cell- 
averaged value of q , and is properly located at the cell centroid; i.e. , at the grid 
point j , k , £ (Ref. 6 ). This geometrical interpretation simplifies the task of devising 
boundary-point difference operators that are consistent with the central differences 
employed at interior points. We consider first the outflow boundary £ = £ max * 

Outflow Boundary. Associated with each point of the outflow boundary is a half-cell, 
as illustrated in Fig. 4-2. The centroid of the half-cell is indicated by the asterisk. 

For the outflow boundary, the only exterior derivative in Eq. (4. 1) is the first term 
in parentheses on the L. H. S. The natural differencing of this term that is compatible 
with the use of (4. 20a) at the adjacent interior point is 


(I % ~ (f j - V 

= 1 


(4.21) 


where the indicated equivalence to the backward difference operator of Eq. (4. 19b) 
follows from Eq. (4. 20b). The time derivative term in (4. 1) must be applied at the 
centroid of the half-cell 



(4.22a) 


4-13 




Fig. 4-2 Sketch of Interior Cells and Boundary Half-Cells Along the 
Coordinate £ 
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where 


= -Ar 
Jkl 


\li. 8 . 

j ] J 
h "i s i [' 


f + w k 6 k 


h - Re 


D - RS 1 *«>] 

1 £<">] - Re 1 ^ k 0 (T? * + | 


(4.17) 


and where it is understood that the spatial derivatives involved in 0^ and co^ are 
differenced as described following Eq. (4.13c), and those involved in 0^, oi^ , R , 
and S are differenced according to Eq. (4.14). 


Each step of the ADI sequence (4. 16) involves the solution of a linear system of equa- 
tions having a block-tridiagonal coefficient matrix. For example, upon applying the 
operator on the L.H.S. of Eq. (4.16a), we obtain 


- Ar 


2 




r! 1 , 2 s j s J - 1 

J J 


(4.18) 


where the subscripts k,l have been suppressed. The block-tridiagonal structure is 
readily apparent, inasmuch as the Jacobian matrices (the blocks) are 5X5 

square matrices. 


However, the system of Eq. (4.18) is incomplete; it involves only J-2 equations among 

the J unknown values of Aq** . The missing equations are those for the boundary points 

j = 1 , J. To close the system, Steger (Ref. 3) and Pulliam and Steger (Ref. 7) simply 

assume Aq|* = Aqjf* = 0 ; i. e. , that the flow conditions at the boundaries do not 

change over a time step. After the sequence (4.16) is completed for interior points, 

— ^ 

they then update the values of q at the boundaries by ad hoc methods that are intended 
to satisfy the boundary conditions at steady state. This time-lagging approach implies 
only a zero-order time accuracy at the boundaries. The approach does not necessarily 
degrade the accuracy of the steady-state solution, depending on the updating procedures, 
but it may retard the rate of convergence to steady state. The numerical experiments 
described in Section 5 suggest that this is the case. 
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Rather than use the time-lagging approach, we follow the approach taken in Ref. 15, 
and formulate implicit difference equations for the boundary grid points that are 
similar to those for interior points, but that embody the boundary conditions of 
Section 2.4. 

4. 2 BOUNDARY POINTS 

In this subsection, we formulate boundary point-difference equations that have the same 
time accuracy as the interior point equations, and that are compatible with the latter 
in their spatial order of accuracy. As indicated in the preceding subsection, the present 
treatment of boundary points differs from that in Refs. 13 , 3 , 6 , and 7 . The present 

treatment for Eqs. (2.20) is patterned after that presented in Ref. 15 for the complete 
Navier-Stokes equations (2.9), in that it provides a fully implicit set of equations for 
the boundary points. 

4.2.1 Time Differencing 


The time differencing for boundary points is the same as that for interior points, as 
given in Eq. (4. 1). However, it is convenient to perform the spatial differencing before 
linearizing the equations in time (cf. subsection 4.1, where the time linearization is 
performed first). With this approach, the boundary conditions enter naturally into the 
difference equations, and the formulation of the implicit operator matrices for the 
viscous terms is effected more easily. 


4. 2. 2 Spatial Differencing 

In addition to the central operators given in Eq. (4.11), we shall employ the classical 

forward difference operator A, and the backward difference operator V. For a mesh 

function f . , the operators for the j direction are defined as follows 
JKl 


AX, = (f lA , - t) 


j jki 


j + 1 j 


(4.19a) 


kf 


V.f., = (f. - f. ) 

j jkf '] j-1 7 


(4. 19b) 


kf 
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The corresponding operators for the k and S. directions are identified by the latter 
indices as subscripts. The subscripts will serve to distinguish the backward difference 
operator from the vector gradient operator (Section 2), both of which are identified by 
the symbol V ; and to distinguish the forward spatial operator from the time operator 
of Eqs. (4.2). 

When Eq. (4. 1) is applied at a point on the boundary of the computational space 4 , V , £ , 
the spatial derivatives in that equation fall into two classes: interior derivatives taken 
in directions that lie in the boundary surface itself, and one-sided exterior derivatives 
taken in the direction normal to the boundary. For example, at the inflow boundary 
4=0 and the outflow boundary 4 = ^ max » the derivatives 9/3 r\ and 3/3 £ are 
interior derivatives along the boundary surface, and 3/34 is an exterior derivative 
in the direction normal to the boundary. All interior derivatives in Eq. (4. 1) can be 
differenced in exactly the same fashion at the boundaries as at interior points of the 
computational space (see subsection 4. 1.2). The exterior derivatives must be treated 
differently. Before doing so, however, it is instructive to interpret geometrically the 
central difference operators that are employed in subsection 4. 1.2. 

As demonstrated in Ref. 6, the use of central difference operators at an interior 
grid point of the computational space is equivalent to applying the finite-volume method 
(Refs. 17, 18 ) to a cubical cell enclosing that grid point. As illustrated in Fig. 4-la 
for a cross-sectional plane £ = const. , half-integer subscripts effectively define the 
faces of a cubical cell of unit volume A4 Arj A£ = 1 that encloses each interior 
point. The dashed lines in the figure delineate the cell boundaries. If Eq. (4. 1) is 
integrated over the volume of cell jkl in the figure and the resulting equation is divided 
by the cell volume A4 At? A£ then, for example, there appears a term that 
represents the first derivative df/d 4 as 


Of784) jk£ 


(f j + l/2,k,f " f j-l/2,k,f )/A4 



A4 = 1 


(4.20a) 
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Fig. 4-la Sketch of Computational Grid in a Plane £ = Constant 
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Fig. 4 -lb Image of Cell j , k , £ in Physical Space 
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In the finite volume method, the value of f at a cell face is defined as the average of 
values at the two neighboring grid points 


f j + l/2,M (f jki + f j + l,M )/2 M j f j + 1/2,M 


(4.20b) 


When the left-hand equality of (4.20b) is substituted into that of (4.20a), the resulting 
equation is identical to Eq. (4. 13a); namely, the central difference that is obtained by 
applying the operators 6^ and in sequence. Similar results hold for all the spatially 
differentiated terms of Eq. (4. 1). Furthermore, by the mean value theorem, the first 
term in the volume integrated version of (4. 1) represents the time derivative of the cell- 

A 

averaged value of q , and is properly located at the cell centroid; i.e. , at the grid 
point j , k , f (Ref. 6 ). This geometrical interpretation simplifies the task of devising 
boundary-point difference operators that are consistent with the central differences 
employed at interior points. We consider first the outflow boundary | = £ max • 


Outflow Boundary. Associated with each point of the outflow boundary is a half-cell, 
as illustrated in Fig. 4-2. The centroid of the half-cell is indicated by the asterisk. 
For the outflow boundary, the only exterior derivative in Eq. (4. 1) is the first term 
in parentheses on the L. H. S. The natural differencing of this term that is compatible 
with the use of (4. 20a) at the adjacent interior point is 


(Hlj ~ < f J - f J-l/2> / < A « /2 > ■ V 

A£ = 1 (4.21) 


where the indicated equivalence to the backward difference operator of Eq. (4. 19b) 
follows from Eq. (4. 20b). The time derivative term in (4. 1) must be applied at the 
centroid of the half-cell 



(4. 22a) 


4-13 



I 


LL'l 

ft 


I * 
l I 

-4- - 

I 

i 


4 = 0 


Fig. 4-2 Sketch of Interior Cells and Boundary Half -Cells Along the 
Coordinate £ 


4-14 



The latter can be related to the corresponding quantities at the neighboring grid points 
within the spatial order of accuracy of (4. 21) by linear interpolation (Refs. 6 and 15) 


Aq* = (3 Aqj + Aqj _ t )/4 

= [I - (1/4) V.] Aqj (4.22b) 

Upon inserting Eqs. (4.21) and (4.22) into Eq. (4. 1) along with the central difference 
operators of subsection 4.1.2 for the interior spatial derivatives 9/9 rj and 9/9£ , 
the time linearization and implicit operator factorization can be performed directly to 
yield an ADI sequence similar to the interior point sequence (4.16) for interior points. 
The first step of the sequence for outflow boundary point is 


[I + ATV. (F (4.23a) 

a = 1/4 At (4.23b) 

where the term on the R.II.S. can be obtained fran the interior point equation (4. 17) 
by the substitution 6. — * in the first term on the right of the latter equation. 

The other steps of the sequence are identical to those given by Eqs. (4. 16b through d). 

The described sequence is valid only for purely supersonic outflow, or for problems 
that involve an external flow. In these cases, the full set of flow equations is used 
at the outflow boundary, and the latter is situated downstream of the nozzle exit 
(see subsection 2.4.6). For purely internal flows, we must modify the ADI sequence 
to account for the boundary condition given in Eq. (2. 31). With the aid of the equation 
of state (2. 6b) and the definition of q given below Eq. (2. 9), the boundary condition 
can be written in the form 

5(q) = 5g - q^/qg + q 3 + q^J - J/y (y - i) = o (4.24) 
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where q. , i = 1,2 5 denote the components of the vector q and J is the 

1 n + 1 

Jacobian of Eq. (2. 10). This boundary condition is to be satisfied at time r 

J n + 1 (q) = 0 (4.25a) 

The latter is linearized in time as follows: 


^n+i = Aq = 0 


We write the linearized boundary condition in the form 


n . a a n 
v Aq - -J 


(4. 25b) 


(4.26a) 


where v denotes the row vector 

d9 

V = = (V 2 /2 , - u , - v, - w, 1) (4. 26b) 

dq 

The linearized algebraic boundary condition is used in place of the u-momentum 
equation, i. e. , in place of the third component of the time-linearized and spatially 
differenced version of Eq. (4. 1). The resulting vector equation can be written in 
the operational form 


|m + Ar [v. (F - cv7) 


W G 




- 


Re 1 


«5 k R 


6 f S ) 


I) A %U 


= r 


Jkf 


(4.27) 


All elements are zero in the third row vector of each of the matrices F, I, G, H and 
of the operator matrices R, S; the matrices otherwise are identical to F^ I, G, H, R, S 
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as defined in subsection 4. 1. 1. The third row vector of the matrix M is v as defined 
in Eq. (4. 26b); the other elements of M are those of the identity matrix I. Finally, 
the third component of the vector on the right of (4.27) is equal to the R. H.S. of 

a 

Eq. (4.26a); the other components are those of in Eq. (4.23a). 

The operator on the left of (4.27) can be factored easily after multiplying the equation 
by M \ and the resulting ADI sequence can be written in the form (Ref. 14 ). 

r*j ' 'r-i 

M + At V (F - al)] Aq Jk£ - r^ (4.28a) 

M + AT (p k <5 k G - Re\R)] Aq Jk£ - M 1 Aq Jkjf (4. 28b) 

[m + At (jn^H - Re^S)] Aq^ - M n Aq^ (4.28c) 

Inflow Boundary. The algebraic inflow boundary conditions on the total pressure, total 

temperature, and transverse velocity components (see subsection 2.4.5) are treated 

in the manner described above for the algebraic outflow boundary condition in 

Eq. (2.31). The implicitly-differenced form of the mass conservation equation is used 

along with the algebraic inflow boundary conditions to form a complete system of 

implicit linearized equations for inflow boundary grid points. The ADI sequence has 

the same form as Eq. (4. 28) and can be obtained from the latter by the substitutions 

J — 1, a — - - a, and V. •— A. . Note that the latter substitution also must be made 

a 1 J 

in the R. H. S. term r . 

Free stream Boundaries. The computational space described in Section 3 is such that 
the freestream boundary conditions of subsection 2 . 4. 3 apply at either of the boundaries 
V ~ 7 7 TV ,_, r » £ - £ » depending upon the nature of the curvilinear coordinate system. 

For the grid geometry depicted in Fig. (3-6a), both are freestream boundaries. We 
describe below the treatment of Eq. (4. 1) for boundary points 77 = r] max (k = K) . 
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The interior derivatives 9/9£ , 9/9? are differenced centrally as for interior and 
points (see subsection 4. 1.2). Backward difference operators are used to represent 
the exterior first derivatives dg/drj and 90^Vcb) in the manner described earlier 
for exterior derivatives at the outflow boundary. The derivative 30' '/dr] is differenced 
as 


90 
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(Arj/2) 


In view of the boundary conditions (2.27b), the spatially differenced form of Eq. (4. 1) 
for the boundary point j, K, £ then becomes 


(I - (1/4) V k > Aq. jKJ + AT (^.f 4 V k g 4 vfp 


n + 1 


AT Re 1 


2 E 


.“V2 f-< r?) 

j,K 


J" 


+ 1 


+ E 


-1 

k 



(4.29) 


where E is the classical shift operator (Ref. 16). For a mesh function f. . , the 

J » k» " 

shift operators in Eq. (4.29) are defined by 
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(4. 30b) 


Upon performing the time-linearization Eqs. (4.3) and (4.5a), the resulting implicit 
difference equation is 
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and 


r jKf 


= - At 
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jK4 


Note that the differential operator d/dt] that enters into the quantities 9 y and R 

[see Eq. (2. 10b) and Eq. (A. 4) of Appendix A] may be replaced by the central difference 

operator 6 , as is done in subsection 4. 1.2 for interior points, because the operator 

- 1/2 k 

shifts the point of application of the operator 6^ backward by half the mesh 
spacing. 


The three algebraic boundary conditions (2. 27a) are used in place of the mass, u- 
momentum, and energy equations, i.e. , in place of the first, second, and fifth 
components of the vector Eq. (4.31). The final ADI sequence for the freestream 
boundary points then can be obtained in a fashion analogous to that employed in 
Eqs. (4.24 ) to (4.28) for outflow boundary points. 

Wall Boundaries. The ADI sequence for wall-boundary points can be obtained in a 
fashion similar to that outlined above for the freestream boundary. As shown in 
Ref. 15 , the algebraic boundary conditions (2.21) replace the three momentum 
equations, namely, the second, third, and fourth components of the vector equation 
(4. 1). When the wall temperature is specified according to Eq. (2.22), the latter 
replaces the energy equation (Ref. 15 ) namely the fifth component of Eq. (4. 1). For 
an adiabatic wall, the boundary conditions (2.21) and (2.23b) simplify the spatially- 
differenced form of the energy equation in much the same way as the freestream 
boundary conditions (2. 27b) contribute to the simplicity of Eqs. (4. 31). 


Symmetry Planes. At a symmetry plane such as the vertical plane rj = 0 for the 
nozzle configuration depicted in Fig. (3. 6a), the exterior derivatives 9/9 r) in 
Eq. (4. 1) are evaluated as described in subsection 2.4.3. In the third component 
of the vector Eq. (4. 1), these derivatives vanish according to Eq. (2.25). In the 
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remaining four components, the differential operator 9/3 17 m 3-y he replaced by the 
forward difference operator A k , which is consistent with Eq. (2.26). 

The treatment of the time derivative term Aq/Ar is different from that used at all 
other types of boundary. At a symmetry boundary, each boundary point is considered 
to lie at the centroid of a full cell enclosing the point. The time -derivative term thus 
is centered at the grid point itself, and the counterpart of Eq. (4,22b) does not apply. 

Once the time-linearization is performed, the resulting implicit operator may be 
factored to obtain the ADI operator sequence corresponding to the interior point 
sequence (4. 16). 

4.3 SMOOTHING 

The algorithm of subsection 4. 1 is neutrally stable (Ref. 13)in that short wavelength 
spatial disturbances are undamped. An explicit smoothing term is appended to the 
R. H. S„ of Eq. (4. 16) to prevent nonlinear instabilities and to provide a dissipative 
mechanism for computing embedded shock waves. Steger (Ref. 3 ) and Pulliam and 
Steger (Ref. 7 ) employ a nonconservative fourth-order smoothing term for each 
coordinate direction. Such higher-order terms do not compromise the second-order 
spatial accuracy of the algorithm. However, the analysis given in Ref. 6 shows that 
these nonconservative smoothing terms generate errors in global conservation that 
degrade the solution. This degradation is confirmed by the numerical experiments 
described in Section 5. This defect is overcome by the conservative smoothing term 
developed by Thomas (Ref. 15). The form of this term for the ij coordinate direction 
is 
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where J is the Jacobian in Eq. (4. 10) and (3 is a constant. The first equation applies 

at those interior points located more than two mesh intervals from the boundaries, 

and the second equation applies at the point j = 2 immediately adjacent to the left 

hand boundary in Fig. 4-2. The corresponding term for the point adjacent to the right 

1 /2 - 1 /2 

hand boundary is obtained from Eq. (4.32b) by the substitution E / - E . . 

No smoothing term is applied at the end points themselves. The smoothing terms in 
Eq. (4. 32) remain valid at boundary points that lie in the surfaces rj = const, or 

S' 

° = const, where the difference operators 6^ represent interior derivatives 3/3£ . 

Similar smoothing terms apply for the other two coordinate directions rj , £ . For grid 
points in the neighborhood of a symmetry plane, the smoothing term for the direction 
normal to the symmetry plane is similar to Eq. (4.32a), with appropriate modification 
to account for the symmetry properties of q. The linear stability boundary of the set 
of smoothing terms is (Kef. 15) 


AT 




l + (Ax./2) max [3(fnJ)/3x.] 

x. 1 

i 


< 1/8 


(4.33) 


where x. , i - 1 , 2 , 3 represents the coordinate directions £ , rj , f , respectively. 

The smoothing coefficients /?. must be selected to satisfy this inequality. 

4.4 PROPERTIES OF THE NUMERICAL METHOD 

A linear stability analysis shows that the difference equations of subsection 4. 1 are 
unconditionally stable for arbitrary values of the time step At (Ref. 19). Since only 
the steady-state solution is of interest in the present application, we point out several 
important features of the method that ensure an accurate steady state solution. 

First, at steady state, the term r on the R. H.S. of Eq. (4. 16a) vanishes [cf. Eq. (2.20)]. 
Because Eqs. (4. 16) yield a direct solution for the change in the flow variable vector 
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over a time step, the fact that r is proportional to At implies that the steady state 
solution is independent of the magnitude of the time step (Ref. 13). The accuracy of 
the solution then depends only on the accuracy with which the spatial derivative terms 
in Eq. (2.20) are represented on the finite-difference grid and, of course, on the 
fineness of the grid. These derivatives are approximated at interior points by 
difference operators whose truncation error is of second order in the spatial grid 
spacing, and can be expected to yield reasonable accuracy for a given grid. The 
difference operators employed at boundary points have a local truncation error that 
is of first order in the grid spacing. However, by examining the entire system of 
difference equations for the computational domain, one can show that global second- 
order accuracy is achieved; i.e. , the global truncation error is of second order 
(Ref. 15). Furthermore, strict global conservation is maintained to within machine 
roundoff error. That is, the difference equations possess the same property of 
global conservation as the partial differential Eq. (2.20). 

The global conservation property of the differential equations spanning the physical flow 
region is that only the fluxes through the boundary of the region contribute to the volume 
integral of the physical variables. For the difference equations given in subsections 
4. 1 and 4. 2, the volume integral of the physical variables is expressed naturally as a 
quadrature sum. When the difference equations are summed with the appropriate 
weights of the quadrature, global conservation is achieved if the physical fluxes make 
a net contribution only at the boundaries, and any ad hoc smoothing terms make no 
contribution. 

It is easy to show that the weighted sum of the difference equations over all interior 
and boundary grid points retains the global conservation property of the differential 
equations when the weights are chosen according to the midpoint quadrature rule. 

The latter is appropriate because it is consistent with the second-order spatial 
accuracy of the central difference operators that are employed at interior points of 
the grid. The smoothing terms presented in subsection 4. 3 are differenced in a 
conservative fashion and obey homogeneous boundary conditions. This ensures that 
the weighted sum of these terms over all grid points forms a telescoping series in 
which the contributions from adjacent interior points cancel, and the contribution from 
each boundary point is zero. Thus, the smoothing terms do not compromise the 
global conservation property of the difference equations. 
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Section 5 

NUMERICAL EXPERIMENTS 


An array of numerical experiments have been conducted to test various aspects of 
the numerical method. The bulk of the experiments have been performed to develop 
and verify methods of computing the boundary conditions that are valid throughout the 
subsonic, transonic, and supersonic flow regimes. Experiments also have been 
performed to explore the effects of explicit smoothing terms on the solution. 

The test problem selected is that of the external laminar flow over a flat plate of 
finite length at a freestream Reynolds number Re^ = 10^ , Prandtl n um ber Pr = 1 , 
and viscosity proportional to temperature, p = T . For the subsonic flow experiments, 
the freestream Mach number is M m = 0.1. The computational grid for this case is 
shown in Figs. 5-1 and 5-2, and features parabolic boundary-layer curvilinear coordi- 
nates. The lateral outer computational boundary is placed approximately ten boundary 
layer thicknesses away from the plate. The grid is stretched exponentially in the z 
direction to resolve the thin viscous layer. The number of grid points in the x , y , 
and z directions is 15 x 5 x 15. The stretching in the z direction places about half 
of the grid points within the boundary layer. 

All test runs employ the implicit adiabatic wall boundary computation technique of 
Ref. 15 (see section 4.2.2). The initial conditions are as follows. Pressure is 
assumed uniform and equal to its freestream value, p = 1 . Velocity is taken from 
the Blasius boundary layer solution for incompressible flow (Ref. 20). Temperature 
is computed from the velocity through the Crocco relation, and density follows from 
the equation of state. 

Each set of experiments is discussed below in a separate subsection. The conclusions 
drawn are summarized separately at the end of each subsection. 
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5.1 OUTFLOW BOUNDARY CONDITION TESTS 


As discussed in section 2.4.5, the quasi-one-dimensional theory of characteristics 
implies that no boundary conditions are to be specified when the x-component of 
outflow velocity is supersonic, u > c . Conditions at the boundary exert no influence 
on the upstream flow, and simple extrapolation methods can be used to compute the 
flow variables at the boundary. An improvement over extrapolation is to compute the 
flow variables from the equations, using "upwind" difference operators for the terms 
that contain spatial derivatives in the direction normal to the boundary (Ref. 6) . Such 
operators are termed "conservative" if the resulting difference equations possess the 
same global conservation property as the differential equations that govern the flow 
(see Section 4.4) . This conservative outflow boundary computation procedure has been 
applied successfully by Thomas and Lombard (Ref. 6) in supersonic external viscous 
flows that contain a subsonic region embedded in the near-wall region of the boundary 
layer. The presence of an outer supersonic region was believed responsible for the 
success of this procedure, inasmuch as characteristics theory suggests that a boundary 
condition is required when the flow is subsonic (see section 2. 4. 6) . 

According to the latter theory, four characteristics are directed from the interior 
toward the boundary, and one emanates from the boundary toward the interior. The 
latter is responsible for the upstream influence of the boundary, and implies that a 
boundary condition is necessary, say, the boundary pressure. The boundary condition, 
together with the four compatibility relations that hold along the other four character- 
istics then could be used to compute the flow variables at points of the boundary where 
u < c . The implicit solution of the compatibility relations in a fashion consistent 
with the implicit method for interior points is much more complicated than the conserva- 
tive outflow boundary computation procedure. Numerical experiments therefore have 
been performed to investigate the validity of the latter procedure for subsonic and 
transonic flows. 

The test case is for the flat plate problem described earlier, with a freestream Mach 
number M =0.1. Freestream boundary conditions p = p = 1 , u = M , v = w = 0 

OO 00 

were imposed at the lateral outer boundary z = z (x) , and held fixed throughout 

max 
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the computation. Conditions at the inflow boundary also were held fixed at their initial 
values. The nonconservative smoothing terms (Refs. 3, 7, and 13) were used. A 
constant time stepsize At = 0.005 was used. This value is twenty times larger than 
the stable limit of explicit numerical schemes as given by the Courant stability 
criterion. That is, the Courant number for the calculation is Co = 20 . 

Figure 5-3 shows the computed surface pressure distribution after 100 steps . The 
pressure peak at the left is probably a result of the simple inflow boundary condition. 
Note the smoothness of the pressure distribution near the trailing edge of the plate 
where the conservative outflow boundary computation procedure is used. This result 
shows that the procedure yields smooth, stable results even when the subsonic region 
is not confined to the near-wall portion of the viscous layer. 

We conclude that the procedure is valid for subsonic and transonic flow, as well as 
for the supersonic external flow problems for which it was employed originally in 
Ref. 6. It is likely that the success of the procedure is due partly to the use of 
conservative boundary-point difference operators, and partly to the fact that the 
viscous flow equations are spatially parabolic in the transverse y , z directions, 
rather than being strictly hyperbolic as assumed in the quasi-one-dimensional theory 
of characteristics. We conjecture that the parabolicity allows information about the 
freestream conditions to be transmitted along the outflow boundary to the wall, and 
obviates the need for the pressure to be specified as a boundary condition . Having 
confirmed its general validity, the conservative outflow boundary computation 
procedure is used in all subsequent numerical computations. 

5.2 FREESTREAM BOUNDARY CONDITIONS 

Although the numerical experiment just described confirmed the validity of the 
conservative outflow boundary condition procedure, it revealed that a truly steady - 
state solution is difficult to obtain for subsonic flow M < 1 if all freestream flow 

CO 

conditions are imposed at the lateral outer boundary z = z . Further computa- 

m 

tion shows that, although the numerical method remains stable for Courant numbers 
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as high as 400, the solution continues to change with time even after many hundreds of 
steps. After many numerical experiments to investigate possible causes of the per- 
sistent unsteadiness in the solution, the evidence became strong that the difficulty lay 
in the use of freestream boundary conditions at the outer boundary. 

It is known from second-order boundary layer theory that the displacement effect of 
the viscous layer near the wall induces a first-order perturbation of the outer inviscid 
flow (Ref. 21) . Specifically, a finite velocity v is induced in the direction normal to 
the wall. Apparently, this induced normal velocity does not decay with distance away 
from the wall in a wholly subsonic flow, and one cannot impose at the outer computa- 
tional boundary the freestream conditions v = w = 0 , no matter how far that boundary 
is placed from the wall. 

A variety of numerical experiments were conducted to develop a valid method for the 
outer boundary. The first of these experimented with a generalization of the conserva- 
tive outflow boundary computation procedure in which the only boundary conditions 
applied were that the viscous stress and heat conduction terms must vanish at the 
boundary because the freestream flow is inviscid. All flow variables then are 
computed implicitly from the mass, momentum, and energy conservation equations. 
However, better results were obtained by imposing freestream boundary conditions 
on the density, streamwise velocity u , and pressure. The transverse velocity 
components v , w are computed from the transverse momentum equations as in the 
generalized outflow boundary computation procedure. 

A steady-state solution was obtained in about 150 steps, using a time stepsize that 
corresponds to a Courant number Co = 200 . The results are displayed in Figs. 5-4 
to 5-11. The first seven of these figures show vertical profiles of the velocity 
components, the density, and the pressure. The velocity profiles at a station midway 
along the plate (Fig. 5-4) are smooth. Note that the vertical velocity v (dashed line) 
remains non-zero all the way to the outer boundary. The density profile of Fig. 5-5 
at the same station is plotted on a greatly expanded scale, because the maximum 
density change between the wall and the freestream boundary is less than 0.2% at the 
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Fig. 5-5 Density Profile at a Station Midway Along Plate, = 0.1 
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Fig. 5-7 Velocity Profiles at Trailing Edge, M =0.1 
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Fig. 5-8 Density Profile at Trailing Edge, M m = 0.1 
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Fig. 5-9 Pressure Profile at Trailing Edge, = 0.1 
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very low freestream Mach number = 0.1 . The pressure is uniform across 

the entire computational region, as shown in Fig. 5-6. 

The profiles are similar at the trailing edge of the plate, which coincides with the 
outflow boundary (Figs. 5-7 to 5-9). Figure 5-10 shows a magnified view of the 
velocity profiles across the near-wall viscous layer at the trailing edge. 

We have found the computed drag coefficient to be a sensitive indicator of convergence 
to steady state. The variation in drag as a function of the number of time steps 
computed is shown in Fig. 5-11. The negative sign is associated with the sign 
convention used to denote the direction of a force component. The steady state is 
seen to be reached in less than 200 time steps, according to the drag computation. 

This was confirmed by the fact that an additional 100 steps produced no further change 
in any flow variable at any grid point, to within machine roundoff error. Furthermore, 
the root mean square residual over all grid points had decayed to the order of the 
machine roundoff error, where the residual is defined as the set of steady-state terms 
in the flow equations. By this definition, the residual should vanish at steady state. 

The asymptotic steady state value of the computed drag coefficient is 
C = 3.97 x io" 3 , about 5 percent lower than the value predicted by the approximate 

Blasius boundary layer solution for incompressible flow (Ref. 21). The Blasius 
value is 


C D = 1.328//Ke^ = 4.20 x 10 -3 (5.1) 

One would expect the same freestream boundary computation procedure to be valid for 
supersonic flow. This was confirmed by a successful computation for a freestream 
Mach number M m = 3 , again at a Reynolds number Re^ = 10 . For a supersonic 

freestream, the viscous interaction induces a shock wave in the predominantly inviscid 
flow outside the viscous wall boundary layer. The computed surface pressure distri- 
bution in Fig. 5-12 shows that the shock-induced elevation in surface pressure is 
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greatest near the leading edge, where the shock is strongest. Note the smooth 
variation in pressure all the way to the trailing edge, where the conservative outflow 
boundary computation procedure is used. 

Figure 5-13 shows the velocity profiles at midplate. The weak embedded shock wave 
is evident in the abrupt decay of the vertical velocity v (dashed line) across the 
shock. The shock is poorly resolved by the exponentially stretched grid used in this 
test case. The shock is again evident in the pressure profile (Fig. 5-14) at the same 
station. The post-shock oscillation is common with such shocks that are "captured” 
by the numerical method, but probably is aggravated by the aforementioned poor 
resolution of the shock jump. 

Similar velocity and pressure profiles at the trailing edge are displayed in Figs. 5-15 
and 5-16. The density profile at the trailing edge is given in Fig. 5-17 to show the 
large density change that occurs across the boundary layer. The density change is a 
consequence of the temperature rise associated with viscous dissipation of kinetic 
energy, which is much greater here than in subsonic flow (cf. Fig. 5-8). Magnified 
views of the velocity, pressure, and density profiles across the inner viscous boundary 
layer region at the trailing edge are given in Figs. 5-18 to 5-20 to show the smoothness 
of the solution, and the good resolution that is achieved with as few as eight grid points 
across the layer. 

The computed drag coefficient history is shown in Fig. 5-21. The Courant number 
based on the time stepsize is about 100 for this case, and steady state is again achieved 
within about 200 time steps. The computed drag coefficient at steady state is 
= 4.54 x 10 -3 . According to laminar boundary-layer theory, the drag is 

independent of Mach number for Pr = 1 and viscosity proportional to temperature, 
and the incompressible Blasius value of Eq. (5.1) applies (Ref. 20). The computed 
drag is about 8 percent larger than the Blasius value. The discrepancy is not 
surprising, in view of the fact that boundary-layer theory neglects the viscous 
interaction, and does not account for the presence of the external shock wave. 
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Fig. 5-13 Velocity Profiles at a Station Midway Along Plate, = 3 
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Fig. 5-14 Pressure Profile at a Station Midway Along Plate, = 3 
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Fig. 5-15 Velocity Profiles at Trailing Edge, = 3 
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Fig. 5-18 Velocity Profiles Across Inner Viscous Layer at Trailing Edge, - 3 
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Fig. 5-20 Density Profile Across Inner Viscous Layer at Trailing Edge, - 3 
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Fig. 5-21 Drag Coefficient as a Function of Number of Time Steps, M m = 3. 
Maximum Courant No. Co = 100 
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From these numerical experiments, we conclude that the described procedures for 
computing the flow at outflow boundaries and at freestream boundaries are valid over 
the whole range of subsonic, transonic, and supersonic flow conditions. The numerical 
experiments played a vital role in the final formulation of the boundary conditions 
presented in Section 2.4, and in the development of the boundary point difference 
equations presented in Section 4.2. 

5.3 EFFECT OF SMOOTHING TERMS ON THE SOLUTION 

Numerical experiments also have been performed to explore the effects of various 
types of explicit smoothing terms on the solution. Both the conservative smoothing 
terms of Section 4.3 and the nonconservative terms of Refs. 3, 7, and 13 have been 
tested. In general, the conservative smoothing terms have been found to be superior 
in two respects. 

First, the nonconservative smoothing terms produce a substantial overshoot in both 
the density and the streamwise velocity component at grid points situated near the 
boundary layer edge, where the curvature in the profiles of these variables is greatest. 
That is, the values exceed the freestream values in this region. For the same value 
of the smoothing coefficient, the conservative smoothing terms produce much less 
overshoot. Some overshoot is still present, as evidenced by the profiles in Figs. 5-4 
to 5-20, which were computed using the conservative smoothing terms. 

Second, the conservative smoothing terms do not cause the difference equations to 
violate the global conservation property of the differential equations that govern the 
flow. This is not true of the nonconservative smoothing terms. The significance 
of the global conservation property became evident when we attempted to compute the 
drag by integrating the momentum defect over all boundaries permeable to the flow. 

The drag computed by this method differed by nearly a factor of two from that obtained 
from a direct integration of the surface stresses because of the fact that the nonconserva- 
tive smoothing terms contribute errors to the global mass, momentum, and energy 
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balances over the flow region. The conservative smoothing terms make no net 
contribution to the global balances. Consequently, the drag computed by momentum 
defect is in agreement with that obtained by integrating surface stresses . 
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Section 6 


TURBULENCE MODELS 

6.1 INTERNAL AND EXTERNAL WALL BOUNDARY LAYERS 

The boundary layer flow along the external and internal nozzle walls, as well as on the 
wedge plug and on the side plates, can be modeled using a two-layer eddy viscosity 
model which has been employed successfully in previous calculations of compressible 
boundary layer flows at subsonic and supersonic Mach numbers. Implicit in this 
model is the assumption of two-dimensional, thin shear layers in local equilibrium; 
i.e. , the turbulent flow in the boundary layer represents local mean flow conditions. 

A relaxation eddy viscosity model is available and will be incorporated into the calcu- 
lations in cases where streamwise pressure gradients lead to flow separation. The 
relaxation model accounts for the lag in turbulence response to rapid changes in mean 
flow conditions. The application of other turbulence models such as those of Refs. 22, 
23, and 24 which entail solving differential equations for the Reynolds stresses, turbu- 
lent kinetic energy, and length scale appear not warranted in the present application. 
Such models do not offer significantly better performance, but only increase computa- 
tional work. 

With these assumptions, the turbulent stresses t in the boundary layer are modeled 
in terms of the eddy viscosity by 


9U 
= dy 


( 6 . 1 ) 


where U is the velocity component parallel to the wall, and y is the coordinate nor- 
mal to and measured from the wall. In the inner region (also known as sublayer), p^ 
is calculated from the Van Driest eddy viscosity formulation with damping in the sub- 
layer which was modified by Cebeci et al. (Ref. 25) for the case of compressible flow 
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with pressure gradient, 
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n t = P k y 


t - 


8U 

9y 


( 6 . 2 ) 


where the damping constant A + depends upon the streamwise pressure gradient as 


A + = 26 [ 1 + 

\ 


\-l/2 
d£\ 
ds J 


(6.3) 
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The constant k ha 


is the friction velocity based upon the wall shear stress and density. 
The constant 1c has the value of 0.40. A slightly different expression for A is ob- 
tained by replacing y in Eq. (6.3) with the sublayer thickness y^ defined by 


+ P y S % 

y s = W 


= 11.8 


(6.4) 


so that 


A + = 26^1+ 11.8 P + ) 


- 1/2 


where P is the pressure gradient parameter 


P + = — 

2 3 ds 
p u r 


(6.5) 


( 6 . 6 ) 


In the outer region of the boundary layer , the eddy viscosity is determined from the 
Clauser defect law 
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(6.7) 
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U c represents the core flow velocity at the edge of the boundary layer. Equation (6.7) 
includes in the denominator Klebanoff's intermittency factor (Ref. 26) in order to ac- 
count for the intermittent nature of turbulent flow in the outer region of the boundary 
layer. The choice of Eqs. (6.2) or (6.7) is made so as to assure continuity of p from 
the inner to the outer region . Equation (6 . 7) involves the core velocity U c and the 
boundary layer thickness 6 , both of which may not be well defined in the nozzle in- 
terior for nonuniform inlet conditions. In this case, a mass-flow averaged core flow 
velocity can be defined to represent the flow velocity outside the boundary layer. The 
boundary layer thickness 6 can be defined as the distance from the wall within which 
the velocity or its gradient dU/dy approaches the corresponding values of the external 
stream to within a specified tolerance. 

6 . 2 INTERNAL CORE FLOW 

The nature of turbulence is characterized by slow response to sudden changes in rates 
of strain, for which reason the turbulence of the nozzle core flow is expected to under- 
go negligible changes during its passage through the nozzle. A constant eddy viscosity 
will be determined for the flow outside the boundary layer from the relationship for a 
round jet according to Schlichting (Ref. 20), 

= x P b U q (6.8) 

where k = 0.0256, b is the nozzle inlet radius, and U Q is the average inlet velocity. 
If the eddy viscosity of the core flow is considerably greater than ^ of the outer 
wall region, it would be advantageous to suppress the intermittency factor in Eq. (6.7) 
since the eddy viscosity in this region is increased by entrainment of turbulence. 
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During the passage of the core flow through the nonaxisymmetric nozzle, lateral redis- 
tribution of temperature by turbulent fluctuations or diffusion will be small, but will be 
taken into account according to specified nonuniformities of the inlet flow. Either tem- 
perature or adiabatic boundary conditions can be specified along the nozzle wall. In 
either case, a constant turbulent Prandtl number of Pr t = 0.9 will be used in the 
energy equation. This value has been shown by Cebeci (Ref. 27) to lead to excellent 
predictions of the boundary layer properties at various Mach numbers up to M - 6. 6. 
The turbulent eddy conductivity k^_ is thus defined by 


Pr t = 


c 

_L_P 


= 0.9 


(6.9) 


The wall boundary layers external to the nozzle belong to computational Region C 
illustrated in Fig. 6-1. They will be computed using the same eddy viscosity models 
described above. Although the sketched geometry does not indicate regions of separated 
flow, exterior flow separation may occur in the region of confluence with the interior 
flow. In such cases, the relaxation eddy viscosity model of Shang and Hankey (Ref. 28) 
can be used to correct the local equilibrium eddy viscosity model for the effects of 
upstream history (nonequilibrium effects) upon the development of turbulence. The 
model computes a nonequilibrium eddy viscosity from local and upstream values, 


M 


n 


+ 0* - M 


H 1 - 


exp 



( 6 . 10 ) 


where the subscript "o" refers to the upstream location where an abrupt change in 
streamwise pressure gradient occurs. The relaxation length A has typically a value 
of 2. 5 to 10. In the calculations of turbulent, separated flow over axisymmetric after- 
body boattail configurations by Holst (Ref. 29), best comparison with experimental 
data was obtained for A = 2. 5. 
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6. 3 FREE SHEAR LAYER AND DEVELOPED JET REGION 

The flow in Region B of Fig. 6-1, outside the boundary layers on the side plates, can 
be characterized as a two-dimensional free shear layer. This flow changes later into 



Fig. 6-1 Definition of Mixing Layer Flow in Region B 

a developed free jet after the mixing layer spreads and envelopes the entire exhaust 
stream. Compressible shear layers of this type with differences in temperature 
and density across the layer have been computed successfully with eddy viscosity models 
where is proportional to the width of the mixing region and to a characteristic 
velocity. As in the case of the internal and external boundary layers, we favor the 
eddy viscosity approach for computational reasons over more complicated models, 
e.g., Varma et al. (Ref. 24). In particular, the model used by Donaldson and Gray 
(Ref. 30) represents a convenient and accurate method for computing the merging of 
the two mixing layers in regions into a single jet, although it may not be necessary to 
carry the computations into the fully developed turbulent region. Ihe eddy viscosity 
of the mixing layers is then given by 

"t = k k < Y - Y i» <U c ■ U ~> < 6 - Ua > 

where U c is the maximum streamwise velocity component of Region B, and Y , Y. 
are coordinates measured from the nozzle centerline as indicated in Fig. 6-1. Using 
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the calculated and experimental data of Ref. 30, k depends upon the Mach number W 
at the point where the shear layer velocity equals (U c + U^j/2, 


0.046 - 0.016 M 

for 

M < 1.3 

0.025 

for 

M > 1.3 


(6. lib) 


The dependence of k upon M implies that the existence of a stationary shock structure 
in a free mixing layer has no first-order effects upon the mixing rate. If the downstream 
computational boundary reaches into the developed jet region, ^ is calculated from 
Eq. (6.11) with Y. = 0 and with U c representing the maximum velocity in the jet. 
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Appendix A 

IMPLICIT OPERATOR MATRICES FOR THE VISCOUS TERMS 


The elements of the operator matrix R in Eq. (4.5a) can be derived most easily as 
follows. We first expand the vector product terms that enter into Eq. (2. 16b) for 9^ 
and eliminate the temperature T by using the equation-of-state (2. 6b). The result is 


9^ = fiJ 


a,u + /3,v + B 0 w 

1 >) K 1 r] 3 r\ 

+ 01 2 V V + 

^3 U 77 + + a 3 W r 1 


« 4 ( e /P) J? + (<* 1 ~a 4 ) (u 2 / 2 )^ + (« 2 ~ o; 4 ) (v 2 /2) r 
+ {ot 3 - a 4 ) (w 2 /2) + ^(uv) + /3 2 (vw) 


(A.l) 


P 3 (wu) 


T? _J 


"l = I Vt) | 2 + *1^/3 a 2 = iVT?l 2 + ^y/3 
= ^xV 3 


/?2 - ^y^ z /3 


Of Q = IVT? I 2 +T] 2 /3 


^3 = Vx /3 


Qf. = (y/Pr)|VT)| 


The corresponding equation for can be obtained from Eqs. A. 1 by making the 

substitutions 6 , r\ — co , f . 


A-l 



Now, R denotes the Jacobian matrix operator 


R 


A 

aq 


where R' 


J 1 R' 

(A. 2a) 

d m 

(A. 2b) 

d~q 


If 9 denotes the i'^ 1 component of the vector 9 ^ and q. denotes the j compo- 
i J 

nent of the vector T, then the elements of the matrix operator R’ are 



If we neglect any dependence of p, 7, or Pr on q, then the matrix elements can be 
computed easily after writing the quantities u,v,w in Eq. (A.l) in terms of the com- 
ponents of 


T ,T 

q* = (P,pu ,pv ,pw ,e) = (q x > q 2 ’ ^3 ’ q 4 ’ q 5> 

u = q 2 /q x v = q g /q 1 w = q 4 /q 1 


(A. 3) 


The operator matrix R' then can be written in the bordered form (cf. Ref. 7) 








o 1 
1 

0 

0 

0 

1 0 

r 1 i 

21 1 




1 0 
1 

1 

r’ 1 

r 31 1 

1 


P 


i o 
1 

t— i 




1 0 
J 

r' 1 

ol 1 

r T 
52 

r ’ 

53 

r' 

54 

1 r' 

1 55 





— 


(A. 4a) 
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where the submatrix operator P may be written as the product of an ordinary matrix 
Q and a scalar differential operator 


p = Q ~Sf( 1 ' / P) (A. 4b) 

where Q is symmetric 

^3 

P 2 (A. 4c) 

“3 

The non-zero border elements of R T are given by 




r M “ “ 4 ^- [<V 2 -e/P)/pj - Oj (u 2 /p) - « 2 ± (v 2 /p) - a, -S. (w 2 /p) 

" 2 “l -p < uv/p > - 2 ' J 2 (vw/p) - 2/3 3 -A <wu/p) 

r 52 ' _r 21 ‘ “4 

r 53 = “ r 31 -“ 4 -^( v /P) 

(A , 4 

r 54 = - r il - “4 ifc < w/p > 
r 55 = “4 W < 1/P > 
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The operator notation employed in Eqs. (A. 4) is to be interpreted as follows. For 
example, the final term that appears in the fifth component of the vector R n Aq of 
Eq. (4.5a) is the term 


( r 5 5 )" ^5 * °4 -k A£ 


(A. 5) 


The latter signifies the differential operator a operating on the product of l/p n 

4 dV 


and AC 


“4 -k < 1/pn > ae $ a di (f) 


(A. 6) 


The elements s„ of the matrix operator S can be obtained from Eqs. (A. 4) by the 
simple substitution 17 — £ . 
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